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

A Düster, H. Bröker, E D Rank

► To cite this version:

A Düster, H. Bröker, E D Rank. The p-version of the finite element method for three-dimensional curved thin walled structures. International Journal for Numerical Methods in Biomedical Engineering, 2001, 10.1002/nme.222. hal-01509450

HAL Id: hal-01509450

https://hal.science/hal-01509450

Submitted on 18 Apr 2017

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished 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.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# The *p*-version of the finite element method for three-dimensional curved thin walled structures

A. Düster†, H. Bröker‡ and E. Rank*,§

Technische Universität München, Fakultät Bauwesen, Lehrstuhl für Bauinformatik, Arcisstr. 21, D-80290 München, Germany

In this paper we present an implementation of a three-dimensional *p*-version for structural problems of solids with almost arbitrarily curved surfaces. Applying the blending function method, complex structures can often be modelled by a few *p*-elements, being the basis for a higher order approximation. Numerical examples will demonstrate, that the *p*-version with anisotropic Ansatz spaces allows to predict the structural behaviour of three-dimensional plates and shells with approximately the same amount of degrees of freedom as in the two-dimensional case, yet significantly more accurate due to the three-dimensional model. Furthermore, it is advantageous to compute complex structures exclusively with three-dimensional discretizations as no special elements are needed to model the transition from dimensionally reduced formulations like plates or shells to fully three-dimensional solid elements. Using the *p*-version with anisotropic Ansatz spaces the whole structure can be efficiently discretized with solid elements, even if the aspect ratio of the elements becomes very large.

KEY WORDS: *p*-version; blending function method; plates; shells

## 1. INTRODUCTION

The *p*-version of the finite element method has turned out to be an efficient discretization strategy for many linear elliptic problems. To name a few consider e.g. the Poisson equation, the Lamé 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–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><img src="image_2.png" /><figcaption>Figure 1. Set of one-dimensional standard and hierarchic shape functions for $p = 1, 2, 3$.</figcaption></figure>

construct (geometrically) refined 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 effects, if the polynomial degree is chosen to be moderately high [6, 9]. This includes, e.g. shear locking effects appearing when Reissner–Mindlin problems are considered or Poisson locking, which plays an important role in elastoplastic problems [10–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 $(\xi, \eta, \zeta)$ as well as for the Ansatz of the three components $\mathbf{u}(\mathbf{x}) = (u_x, u_y, u_z)^T$ of the displacement field. The use of anisotropic Ansatz spaces leads to very efficient 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 significantly 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 efficiency of the *p*-version for thin walled structures will then be demonstrated by several numerical examples in Section 4.

## 2. THE *p*-VERSION OF THE FINITE ELEMENT METHOD

### 2.1. *The one-dimensional hierarchic basis*

Following Szabó and Babuška [7] it will first be shown how hierarchic basic functions can be implemented up to any desired polynomial degree. Let us start with the *standard* finite element basis (nodal basis) in one dimension on a *standard element* $\Omega_{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 difference
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 one-dimensional hierarchic shape functions

$$N_1(\xi) = \frac{1}{2}(1 - \xi) \qquad (1)$$

$$N_2(\xi) = \frac{1}{2}(1 + \xi) \qquad (2)$$

$$N_i(\xi) = \phi_{i-1}(\xi), \quad i = 3, 4, \dots, p+1 \qquad (3)$$

with

$$\phi_j(\xi) = \sqrt{\frac{2j-1}{2}} \int_{-1}^{\xi} P_{j-1}(t) dt = \frac{1}{\sqrt{4j-2}} (P_j(\xi) - P_{j-2}(\xi)), \quad j=2,3,\dots \qquad (4)$$

where $P_j(\xi)$ are the well-known Legendre polynomials

$$P_k(\xi) = \frac{1}{2^k k!} \frac{d^k}{d\xi^k} (\xi^2 - 1)^k, \quad k = 0, 1, \dots \qquad (5)$$

The linear functions $N_1(\xi), N_2(\xi)$ are called nodal shape functions or nodal modes. Because

$$N_i(-1) = N_i(1) = 0, \quad i = 3, 4, \dots \qquad (6)$$

the functions $N_i(\xi)$, $i=3,4,\dots$ are called internal shape functions, internal modes or bubble modes. The orthogonality property of Legendre polynomials implies

$$\int_{-1}^{1} \frac{dN_i}{d\xi} \frac{dN_j}{d\xi} d\xi = \delta_{ij}, \quad i,j \ge 3 \qquad (7)$$

In Reference [18] it was shown, that the condition number of the stiffness 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 stiffness 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, stiffness matrices corresponding to polynomial order 1 to $p-1$ are submatrices of the stiffness 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.

## 2.2. 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ó and Babuška [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.
<figure><img src="image_3.png" /><figcaption>Figure 2. Standard hexahedral element $\Omega_{st}^h$: definition of nodes, edges, faces and polynomial degree.</figcaption></figure>

* Hexahedral elements are especially well suited for thin walled structures. One local variable can be identified to correspond with the thickness direction. Therefore it is possible to choose the polynomial degree in thickness direction different from those in in-plane direction.

* The numerical integration of hexahedral elements can be easily performed using a Gaussian quadrature scheme.

The three-dimensional shape functions can be classified into four groups:

1. *Nodal modes:*

$$N_{1,1,1}^{N_i}(\xi, \eta, \zeta) = \frac{1}{8}(1 + \xi\eta)(1 + \eta\zeta)(1 + \zeta\xi), \quad i = 1, \dots, 8.$$

The nodal modes are the standard trilinear shape functions, well known from the isoparametric eight-noded brick element. $(\xi_i, \eta_i, \zeta_i)$ are the local co-ordinates of the ith node of the standard element (see Figure 2).

2. *Edge modes:*

These modes are defined for each individual edge separately. If we consider e.g. edge 1, the corresponding edge modes read: $N_{i,1,1}^{E1}(\xi, \eta, \zeta) = \frac{1}{4}(1-\eta)(1-\zeta)\phi_i(\xi)$.

3. *Face modes:*

These modes are defined for each individual face separately. If we consider e.g. face 1, the corresponding face modes read: $N_{i,j,1}^{F1}(\xi, \eta, \zeta) = \frac{1}{2}(1-\zeta)\phi_i(\xi)\phi_j(\eta)$.

4. *Internal modes:*

$$N_{i,j,k}^{\text{int}}(\xi, \eta, \zeta) = \phi_i(\xi)\phi_j(\eta)\phi_k(\zeta).$$

The internal modes are purely local and vanish at the faces of the hexahedral element.

The indices *i*, *j*, *k* of the shape functions denote the polynomial degrees in the local directions $\xi, \eta, \zeta$.

Three different types of Ansatz spaces have been implemented: the *trunk space* $S_{\text{ts}}^{p_{\xi}, p_{\eta}, p_{\zeta}}(\Omega_{\text{st}}^{h})$, the *tensor product space* $S_{\text{ps}}^{p_{\xi}, p_{\eta}, p_{\zeta}}(\Omega_{\text{st}}^{h})$ and the *anisotropic tensor product space* $S^{p, p, q}(\Omega_{\text{st}}^{h})$. A detailed description of the three Ansatz spaces is given in Appendix A. For the definition of the spaces $S_{\text{ts}}^{p_{\xi}, p_{\eta}, p_{\zeta}}(\Omega_{\text{st}}^{h})$ and $S^{p, p, q}(\Omega_{\text{st}}^{h})$ see also Szabó and Babuška [7].
The polynomial degree for the Ansatz spaces $\mathcal{S}_{\text{ts}}^{p_{\xi}, p_{\eta}, p_{\zeta}}(\Omega_{\text{st}}^{\text{h}})$ and $\mathcal{S}_{\text{ps}}^{p_{\xi}, p_{\eta}, p_{\zeta}}(\Omega_{\text{st}}^{\text{h}})$ can be varied in each local direction separately (see Figure 2). The difference between the trunk space and the tensor product space is relevant for the face modes and the internal modes. For explanation we first consider the face modes, e.g. the modes for face 1. Indices $i, j$ denote the polynomial degrees of the face modes in $\xi$- and $\eta$-direction, respectively.

$$ \text{Face modes (face 1): } N_{i,j,1}^{F1}(\xi, \eta, \zeta) = \frac{1}{2}(1 - \zeta)\phi_i(\xi)\phi_j(\eta) $$

$$ \begin{array}{l@{\qquad}l} \text{Trunk space} & \text{Tensor product space} \\ i = 2, \dots, p_{\xi} - 2 & i = 2, \dots, p_{\xi} \\ j = 2, \dots, p_{\eta} - 2 & j = 2, \dots, p_{\eta} \\ i + j = 4, \dots, \max\{p_{\xi}, p_{\eta}\} & \end{array} $$

The definition of the set of internal modes is very similar. Indices $i, j, k$ denote now the polynomial degrees in the three local directions $\xi, \eta, \zeta$.

$$ \text{Internal modes: } N_{i,j,k}^{\text{int}}(\xi, \eta, \zeta) = \phi_i(\xi)\phi_j(\eta)\phi_k(\zeta) $$

$$ \begin{array}{l@{\qquad}l} \text{Trunk space} & \text{Tensor product space} \\ i = 2, \dots, p_{\xi} - 4 & i = 2, \dots, p_{\xi} \\ j = 2, \dots, p_{\eta} - 4 & j = 2, \dots, p_{\eta} \\ k = 2, \dots, p_{\zeta} - 4 & k = 2, \dots, p_{\zeta} \\ i + j + k = 6, \dots, \max\{p_{\xi}, p_{\eta}, p_{\zeta}\} & \end{array} $$

The Ansatz space $\mathcal{S}^{p, p, q}(\Omega_{\text{st}}^{\text{h}})$ defines 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 $\xi$ and $\eta$-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 $\mathcal{S}_{\text{ts}}^{p_{\xi}, p_{\eta}, p_{\zeta}}(\Omega_{\text{st}}^{\text{h}})$ with $p = p_{\xi} = p_{\eta}$. $q$ defines the degree of all shape functions being of higher order in $\zeta$-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 $\mathcal{S}_{\text{ps}}^{p_{\xi}, p_{\eta}, p_{\zeta}}(\Omega_{\text{st}}^{\text{h}})$ with $p = p_{\xi} = p_{\eta}$ and $q = p_{\zeta}$. Considering a polynomial degree $p = q = p_{\xi} = p_{\eta} = p_{\zeta}$ one observes, that the number of internal modes of $\mathcal{S}^{p, p, q}(\Omega_{\text{st}}^{\text{h}})$ is higher than the one of the trunk space $\mathcal{S}_{\text{ts}}^{p_{\xi}, p_{\eta}, p_{\zeta}}(\Omega_{\text{st}}^{\text{h}})$ and less than the one of the tensor product space $\mathcal{S}_{\text{ps}}^{p_{\xi}, p_{\eta}, p_{\zeta}}(\Omega_{\text{st}}^{\text{h}})$ (see Appendix A). Furthermore, it can be noted, that the internal modes of the space $\mathcal{S}^{p, p, q}(\Omega_{\text{st}}^{\text{h}})$ are polynomials in which the direction $\zeta$ is preferred.

Due to the built-in anisotropic behaviour of the Ansatz space $\mathcal{S}^{p, p, q}(\Omega_{\text{st}}^{\text{h}})$ 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 $\zeta$ 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 differently. The numerical examples in Section 4 will demonstrate that anisotropic Ansatz spaces lead to efficient discretizations.
<figure><img src="image_4.png" /><figcaption>Figure 3. Modelling plate-like structures with hexahedral elements.</figcaption></figure>

Our implementation of the p-version allows not only to vary the polynomial degree for the three different local directions but also to choose a different degree for each primary variable. The following two examples illustrate how to define a polynomial degree template p for a structural problem with three primary variables $u = (u_x, u_y, u_z)^T$:

* Let

$$ \mathbf{p} = \begin{pmatrix} \xi \\ \eta \\ \zeta \end{pmatrix} \begin{pmatrix} u_x & u_y & u_z \\ 1 & 4 & 7 \\ 2 & 5 & 8 \\ 3 & 6 & 9 \end{pmatrix} $$

Considering the trunk space $\mathcal{S}_{\text{ts}}^{p_\xi, p_\eta, p_\zeta}(\Omega_{\text{st}}^h)$ the polynomial degree template defines the Ansatz for the displacement field

$$ u_x \in \mathcal{S}_{\text{ts}}^{1,2,3}(\Omega_{\text{st}}^h), \quad u_y \in \mathcal{S}_{\text{ts}}^{4,5,6}(\Omega_{\text{st}}^h) \quad \text{and} \quad u_z \in \mathcal{S}_{\text{ts}}^{7,8,9}(\Omega_{\text{st}}^h) $$

In case of the tensor product space $\mathcal{S}_{\text{ps}}^{p_\xi, p_\eta, p_\zeta}(\Omega_{\text{st}}^h)$ the definition is to be understood analogously.

* Setting

$$ \mathbf{p} = \begin{pmatrix} \xi \\ \eta \\ \zeta \end{pmatrix} \begin{pmatrix} u_x & u_y & u_z \\ 1 & 3 & 5 \\ 1 & 3 & 5 \\ 2 & 4 & 6 \end{pmatrix} $$

will lead in case of the Ansatz space $\mathcal{S}^{p, p, q}(\Omega_{\text{st}}^h)$ to the Ansatz:

$$ u_x \in \mathcal{S}^{1,1,2}(\Omega_{\text{st}}^h), \quad u_y \in \mathcal{S}^{3,3,4}(\Omega_{\text{st}}^h) \quad \text{and} \quad u_z \in \mathcal{S}^{5,5,6}(\Omega_{\text{st}}^h) $$

This polynomial degree template p can be defined for each individual element or for a group of elements. Our finite element code automatically takes care of interelement continuity of the Ansatz. If adjacent elements have different definitions of p, then always the highest degree for the common edge and face modes is chosen.
<figure><img src="image_5.png" /><figcaption>Figure 4. Number of degrees of freedom for 3D Lamé problems for one hexahedral element and uniform p.</figcaption></figure>

In Figure 4 (left part) the number of degrees of freedom of one hexahedral element with isotropic polynomial degree, i.e. for

$$ \mathbf{p} = \begin{pmatrix} \xi & u_x & u_y & u_z \\ \eta & p & p & p \\ \zeta & p & p & p \end{pmatrix} $$

for a three-dimensional Lamé 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. $\mathcal{S}_{\text{ps}}^{P\zeta, P\eta, P\zeta}(\Omega_{\text{st}}^{\text{h}})$ 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 stiffness matrix is strongly reduced. Several authors [19–21] have investigated these observations in detail, interpreting the internal mode condensation as a preconditioning procedure for iterative solvers.

### 3. THE BLENDING FUNCTION METHOD

One important difference between h- and p-version finite 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 finite 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 finite element analysis and taking advantage of the blending function method, some basic concepts of element matrix
<figure><img src="image_6.png" /><figcaption>Figure 5. Mapping by the blending function method.</figcaption></figure>

computation should be reviewed in the light of this mapping technique. The definition of an element stiffness matrix and an element volume load vector

$$ \mathbf{K}^{e} = \int_{-1}^{1} \int_{-1}^{1} \int_{-1}^{1} \mathbf{M}_{\mathrm{K}}(\mathbf{N}, \mathbf{C}, \mathbf{Q}, \mathbf{J}^{-1}, |\mathbf{J}|) \mathrm{d}\xi \mathrm{d}\eta \mathrm{d}\zeta \quad (8) $$

$$ \mathbf{b}^{e} = \int_{-1}^{1} \int_{-1}^{1} \int_{-1}^{1} \mathbf{M}_{\mathrm{b}}(\mathbf{N}, \mathbf{b}, \mathbf{Q}, \mathbf{J}^{-1}, |\mathbf{J}|) \mathrm{d}\xi \mathrm{d}\eta \mathrm{d}\zeta \quad (9) $$

is given by an integral of matrix functions $\mathbf{M}_{\mathrm{K}}$, $\mathbf{M}_{\mathrm{b}}$, depending on:

* shape functions, $\mathbf{N}$,

* material properties, $\mathbf{C}$,

* body force, $\mathbf{b}$,

* mapping function, $\mathbf{Q}$,

* Jacobian matrix, with its inverse and determinant, $\mathbf{J}, \mathbf{J}^{-1}, |\mathbf{J}|$.

Let now $\mathbf{M}$ be a generic matrix function corresponding to $\mathbf{M}_{\mathrm{K}}$ or $\mathbf{M}_{\mathrm{b}}$. In general integration (8) and (9) is performed numerically in the local co-ordinate system $(\xi, \eta, \zeta)$ of a standard hexahedral element domain $\Omega_{\mathrm{st}}^{\mathrm{h}} = [(-1, 1) \times (-1, 1) \times (-1, 1)]$. Using e.g. a Gaussian quadrature the integral is replaced by a weighted sum

$$ \int_{-1}^{1} \int_{-1}^{1} \int_{-1}^{1} \mathbf{M}(\mathbf{N}, \mathbf{C}, \mathbf{b}, \mathbf{Q}, \mathbf{J}^{-1}, |\mathbf{J}|) \mathrm{d}\xi \mathrm{d}\eta \mathrm{d}\zeta \\ \approx \sum_{l} \sum_{m} \sum_{n} \mathbf{M}(\mathbf{N}, \mathbf{C}, \mathbf{b}, \mathbf{Q}, \mathbf{J}^{-1}, |\mathbf{J}|) |_{(\xi_l, \eta_m, \zeta_n)} \mathbf{W}_{\xi_l} \mathbf{W}_{\eta_m} \mathbf{W}_{\zeta_n} \quad (10) $$

Therefore, the terms $\mathbf{N}, \mathbf{C}, \mathbf{b}, \mathbf{Q}, \mathbf{J}^{-1}, |\mathbf{J}|$ have to be computed only at integration points $(\xi_l, \eta_m, \zeta_n)$. In the following, we will consider the blending function method and show, that $\mathbf{Q}$ and $\mathbf{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ályfalvi and Szabó [15]. Consider a hexahedral element, as pictured in Figure 5. $\mathbf{X}_i = (X_i, Y_i, Z_i)$, $i = 1, \dots, 8$ denote the global co-ordinates of the nodes. $\mathbf{E}_i = (E_{ix}, E_{iy}, E_{iz})$, $i = 1, \dots, 12$ are functions which depend on local co-ordinates $(\xi, \eta, \zeta)$ and describe the shape of each edge. $\mathbf{F}_i = (F_{ix}, F_{iy}, F_{iz})$, $i = 1, \dots, 6$ denote the functions
describing the shape of each face. The mapping function $\mathbf{Q}^e(\xi, \eta, \zeta)$ from local $\xi = (\xi, \eta, \zeta)^T$ to global co-ordinates $\mathbf{x} = (x, y, z)^T$ is obtained by

$$ \mathbf{x} = \mathbf{Q}^e(\xi, \eta, \zeta) = \sum_{i=1}^{8} N_{1,1,1}^{N_i}(\xi, \eta, \zeta)\mathbf{X}_i + \sum_{i=1}^{6} \mathbf{f}_i(\xi, \eta, \zeta) - \sum_{i=1}^{12} \mathbf{e}_i(\xi, \eta, \zeta) \quad (11) $$

The first 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 $\mathbf{f}_6(\xi, \eta)$ describes the parametric mapping of the local $(\xi, \eta)$-plane to the surface of the element:

$$ \mathbf{f}_6(\xi, \eta, \zeta) := \left( \mathbf{F}_6(\xi, \eta) - \frac{1}{4}((1-\xi)(1-\eta)\mathbf{X}_5 + (1+\xi)(1-\eta)\mathbf{X}_6 + (1+\xi)(1+\eta)\mathbf{X}_7 + (1-\xi)(1+\eta)\mathbf{X}_8) \right) \left( \frac{1+\zeta}{2} \right) \quad (12) $$

At face 6, where $\zeta=1$ the blending term $(1+\zeta)/2$ equals 1 and therefore $\mathbf{f}_6(\xi, \eta, \zeta=1)$ describes the difference between the curved and the bilinear face co-ordinates $(\xi, \eta)$. Due to the blending term it is guaranteed, on the other hand, that this difference (i.e. the function $\mathbf{f}_6(\xi, \eta, \zeta)$) decreases linearly to the opposite face 1, where $\zeta = -1$ such that $\mathbf{f}_6(\xi, \eta, -1) = 0$.

The third term in Equation (11) corresponds to the edge blending. Considering e.g. edge 1 we have

$$ \mathbf{e}_1(\xi, \eta, \zeta) := \left( \mathbf{E}_1(\xi) - \frac{(1-\xi)\mathbf{X}_1 + (1+\xi)\mathbf{X}_2}{2} \right) \left( \frac{1-\eta}{2} \right) \left( \frac{1-\zeta}{2} \right) \quad (13) $$

$\mathbf{e}_i(\xi, \eta, \zeta)$, $i=1, \dots, 12$ denote the difference 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 $\mathbf{f}_i(\xi, \eta, \zeta)$, $i=1, \dots, 6$ and $\mathbf{e}_i(\xi, \eta, \zeta)$, $i=1, \dots, 12$ in Equation (11) by inserting Equations (B2), (B1) of Appendix B and rearranging terms we finally obtain

$$ \begin{aligned} \mathbf{x} &= \mathbf{Q}^e(\xi, \eta, \zeta) \\ &= \frac{1}{2}((1-\zeta)\mathbf{F}_1(\xi, \eta) + (1-\eta)\mathbf{F}_2(\xi, \zeta) + (1+\zeta)\mathbf{F}_3(\eta, \zeta) \\ &\quad + (1+\eta)\mathbf{F}_4(\xi, \zeta) + (1-\zeta)\mathbf{F}_5(\eta, \zeta) + (1+\zeta)\mathbf{F}_6(\xi, \eta)) \\ &\quad - \frac{1}{4}((1-\zeta)(1-\eta)\mathbf{E}_1(\xi) + (1-\zeta)(1+\zeta)\mathbf{E}_2(\eta) + (1-\zeta)(1+\eta)\mathbf{E}_3(\zeta) \\ &\quad + (1-\zeta)(1-\zeta)\mathbf{E}_4(\eta) + (1-\zeta)(1-\eta)\mathbf{E}_5(\zeta) + (1+\zeta)(1-\eta)\mathbf{E}_6(\zeta) \\ &\quad + (1+\zeta)(1+\eta)\mathbf{E}_7(\zeta) + (1-\zeta)(1+\eta)\mathbf{E}_8(\zeta) + (1+\zeta)(1-\eta)\mathbf{E}_9(\zeta) \end{aligned} $$
$$ 
\begin{aligned}
&+ (1 + \zeta)(1 + \xi)\mathbf{E}_{10}(\eta) + (1 + \zeta)(1 + \eta)\mathbf{E}_{11}(\xi) + (1 + \zeta)(1 - \xi)\mathbf{E}_{12}(\eta) \\
&+ N_{1,1,1}^{N_1}(\xi, \eta, \zeta)\mathbf{X}_1 + N_{1,1,1}^{N_2}(\xi, \eta, \zeta)\mathbf{X}_2 + N_{1,1,1}^{N_3}(\xi, \eta, \zeta)\mathbf{X}_3 \\
&+ N_{1,1,1}^{N_4}(\xi, \eta, \zeta)\mathbf{X}_4 + N_{1,1,1}^{N_5}(\xi, \eta, \zeta)\mathbf{X}_5 + N_{1,1,1}^{N_6}(\xi, \eta, \zeta)\mathbf{X}_6 \\
&+ N_{1,1,1}^{N_7}(\xi, \eta, \zeta)\mathbf{X}_7 + N_{1,1,1}^{N_8}(\xi, \eta, \zeta)\mathbf{X}_8
\end{aligned} 
\quad (14) $$

Consider now the Jacobian matrix

$$ \mathbf{J}(\xi, \eta, \zeta) = \begin{pmatrix} \frac{\partial x}{\partial \xi} & \frac{\partial y}{\partial \xi} & \frac{\partial z}{\partial \xi} \\ \frac{\partial x}{\partial \eta} & \frac{\partial y}{\partial \eta} & \frac{\partial z}{\partial \eta} \\ \frac{\partial x}{\partial \zeta} & \frac{\partial y}{\partial \zeta} & \frac{\partial z}{\partial \zeta} \end{pmatrix} \quad (15) $$

It contains the derivatives of the mapping function $\mathbf{Q}^e$ with respect to the local co-ordinates $\xi$, $\eta$ and $\zeta$. From (14) it can be readily seen, that all coefficients of $\mathbf{J}$ depend linearly on the following three groups of geometric information:

* the co-ordinates of the nodes $\mathbf{X}_i$ for $i = 1, \dots, 8$,

* the 12 tangent vectors

$$ 
\begin{aligned}
\frac{\partial \mathbf{E}_i}{\partial r} \quad &\text{where } r = \xi \quad \text{for } i = 1, 3, 9, 11 \\
&\quad r = \eta \quad \text{for } i = 2, 4, 10, 12 \\
&\quad r = \zeta \quad \text{for } i = 5, 6, 7, 8
\end{aligned} 
\quad (16) $$

* the six tangential planes

$$ 
\begin{aligned}
\left( \frac{\partial \mathbf{F}_i}{\partial r}, \frac{\partial \mathbf{F}_i}{\partial s} \right) \quad &\text{where } r = \xi, s = \eta \quad \text{for } i = 1, 6 \\
&\quad r = \xi, s = \zeta \quad \text{for } i = 2, 4 \\
&\quad r = \eta, s = \zeta \quad \text{for } i = 3, 5
\end{aligned} 
\quad (17) $$

Therefore, the Jacobian matrix at any interior point $\mathbf{x}(\xi, \eta, \zeta)$ of the element can be computed only from nodal, edge and surface data and from the local co-ordinates $(\xi, \eta, \zeta)$ themselves. For given $(\xi, \eta, \zeta)$, 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 $\mathbf{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 $\mathbf{F}_i$ and $\mathbf{E}_i$ are defined as
<figure><img src="image_7.png" /><figcaption>Figure 6. Geometric information of edges and faces with influence on the Jacobian at an integration point.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Figure 7. Parameter values of the objects.</figcaption></figure>

composed mappings in general. Figure 7 sketches these mappings for the upper face $\mathbf{F}_6$ of the hexahedron (see also Figure 5). Defining $\mathbf{u} := (u, v)^T = (u(\xi, \eta, \zeta = 1), v(\xi, \eta, \zeta = 1))^T$ and $\mathbf{x} := (x, y, z)^T = (x(u, v), y(u, v), z(u, v))^T \mathbf{F}_6$ is given by

$$ \mathbf{F}_6 = \mathbf{x}(\mathbf{u}(\xi, \eta, \zeta = 1)) \quad (18) $$

Assuming that the mapping $\mathbf{x}$ is bijective, we denote by $\mathbf{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 $\mathbf{x}^{-1}$. As in general edges $E_i$ of a hexahedron are obtained by arbitrary geometric operations, their image $\mathbf{x}^{-1}(E_i)$ is not available in closed form and must be approximated by a pointwise
<figure><img src="image_9.png" /><figcaption>Figure 8. Client-Server structure for computing element stiffness matrices and load vectors.</figcaption></figure>

iterative backtransformation. Once this approximation is constructed, the mapping **u** from $(\xi, \eta)$- to $(u, v)$-plane can be defined 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 definition of shape functions is possible. Figure 8 sketches a client-server structure of this computation. At each integration point $(\xi_l^e, \eta_m^e, \zeta_n^e)$ the mapped co-ordinates $\mathbf{Q}^e(\xi_l^e, \eta_m^e, \zeta_n^e)$ and the Jacobian matrix $\mathbf{J}(\xi_l^e, \eta_m^e, \zeta_n^e)$ 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 efficiency we do not invoke 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 effort due to the computation of the mapping function $\mathbf{Q}^e$, the Jacobian matrices $\mathbf{J}$ and the
two-dimensional discretization

<figure><img src="image_10.png" /></figure>

p-version
with 25
quadrilaterals

three-dimensional discretizations

<figure><img src="image_11.png" /></figure>

p-version
with 25
hexahedrals

<figure><img src="image_12.png" /></figure>

h-version
with up to 22496
hexahedrals

<caption>Figure 9. Two- and three-dimensional discretizations of a clamped plate under constant load.</caption>

interprocess communication is very small. Considering e.g. the numerical example presented in Section 4.4, the additional time is only about 2 per cent.

Summarizing, this blending function technique for mapping the geometry of *p*-elements offers the possibility to completely separate all geometric computations involved in a finite 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 different computer, is directly linked to a finite 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 offers the advantage of using all state-of-the-art CAD-techniques like geometric editing or parametric design in a finite element analysis, immediately. The increase of efficiency 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.

# 4. NUMERICAL EXAMPLES

## 4.1. Clamped plate under uniform load

In this first 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* definition 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 = 3000000$ being Young's modulus and $\nu = 0.2$ denoting Poisson's ratio. Considering the two-dimensional discretization based on the Reissner-Mindlin plate theory the warping coefficient accounting for non-uniform shear distributions is chosen to be $\kappa = 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
<figure><img src="image_13.png" /><figcaption>Figure 10. Comparison of the two- and three-dimensional solution for a clamped plate under constant load.</figcaption></figure>

quadrilateral *p*-elements utilizing the tensor product space $\mathcal{S}_{\text{ps}}^{p,p}(\Omega_{\text{st}}^{q})$, 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 refined 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 $\mathcal{S}_{\text{ts}}^{p,p,p}(\Omega_{\text{st}}^{h})$ 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 difference 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 influence of the numerical error is negligible. The stress components $\sigma_{xx}, \sigma_{yy}$ and $\sigma_{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 $\sigma_{zz}$ and $\sigma_{yz}, \sigma_{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 effects, 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 $\mathcal{S}_{\text{ts}}^{p_{\zeta}, p_{\eta}, p_{\zeta}}(\Omega_{\text{st}}^{h})$ ($p = p_{\zeta} = p_{\eta} = p_{\zeta}$ is chosen in this example uniformly for all components) supplies with $p=4$ and
<figure><img src="image_14.png" /><figcaption>Figure 11. Comparison of the h- and p-version for a clamped plate under constant load.</figcaption></figure>

1230 degrees of freedom an accuracy which is not reached by the h-version, even with 98913 unknowns.

Furthermore it is obvious, that significant effort can be saved, if anisotropic Ansatz spaces are used in the *p*-version. To get an indication for an economic *p*-distribution for plate-like structures, we will first investigate the dependence of strain energy upon the polynomial degree template **p** using all three Ansatz spaces $\mathcal{S}^{p, p,q}(\Omega_{\text{st}}^{\text{h}})$, $\mathcal{S}_{\text{ts}}^{p_{\zeta}, p_{\eta}, p_{\zeta}}(\Omega_{\text{st}}^{\text{h}})$, $\mathcal{S}_{\text{ps}}^{p_{\zeta}, p_{\eta}, p_{\zeta}}(\Omega_{\text{st}}^{\text{h}})$. 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.

$$ \begin{array}{l} \xi \begin{pmatrix} u_x & u_y & u_z \\ p+1 & p+1 & p \\ p+1 & p+1 & p \\ p+1 & p+1 & p \end{pmatrix} \\ \eta \begin{pmatrix} u_x & u_y & u_z \\ p & p & p+1 \\ p & p & p+1 \\ p & p & p+1 \end{pmatrix} \\ \zeta \begin{pmatrix} u_x & u_y & u_z \\ p & p & p+1 \\ p & p & p+1 \\ p & p & p+1 \end{pmatrix} \end{array} \quad \text{vs} \quad \begin{array}{l} \xi \begin{pmatrix} u_x & u_y & u_z \\ p & p & p+1 \\ p & p & p+1 \\ p & p & p+1 \end{pmatrix} \\ \eta \begin{pmatrix} u_x & u_y & u_z \\ p & p & p+1 \\ p & p & p+1 \\ p & p & p+1 \end{pmatrix} \\ \zeta \begin{pmatrix} u_x & u_y & u_z \\ p & p & p+1 \\ p & p & p+1 \\ p & p & p+1 \end{pmatrix} \end{array} $$

We will also investigate the difference between a raise of the polynomial degree in in-plane $(\xi, \eta)$ and in transversal $(\zeta)$ direction, i.e.

$$ \begin{array}{l} \xi \begin{pmatrix} u_x & u_y & u_z \\ p+1 & p+1 & p+1 \\ p+1 & p+1 & p+1 \\ p & p & p \end{pmatrix} \\ \eta \begin{pmatrix} u_x & u_y & u_z \\ p & p & p+1 \\ p & p & p+1 \end{pmatrix} \\ \zeta \begin{pmatrix} u_x & u_y & u_z \\ p & p & p+1 \\ p & p & p+1 \end{pmatrix} \end{array} \quad \text{vs} \quad \begin{array}{l} \xi \begin{pmatrix} u_x & u_y & u_z \\ p & p & p \\ p & p & p \\ p+1 & p+1 & p+1 \end{pmatrix} \\ \eta \begin{pmatrix} u_x & u_y & u_z \\ p & p & p \\ p & p & p \\ p+1 & p+1 & p+1 \end{pmatrix} \\ \zeta \begin{pmatrix} u_x & u_y & u_z \\ p & p & p \\ p & p & p \\ p+1 & p+1 & p+1 \end{pmatrix} \end{array} $$

From the upper parts of Figures 12–14, it is evident, that a discretization supplies less error, if the polynomial degree is raised rather for the deflection $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 deflection $u_z$. Therefore, an accurate
<figure><img src="image_15.png" /><figcaption>Figure 12. Ansatz space $\mathcal{S}_{\text{ts}}^{p_{\xi}, p_{\eta}, p_{\zeta}(\Omega_{\text{st}}^{h})}$ with $p = 1, \dots, 8$.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Figure 13. Ansatz space $\mathcal{S}_{\text{ps}}^{p_{\xi}, p_{\eta}, p_{\zeta}(\Omega_{\text{st}}^{h})}$ with $p = 1, \dots, 6$.</figcaption></figure>

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 $\xi, \eta$ or $\zeta$-direction, it has to be mentioned, that the p-version mesh is constructed such that the local coordinates $\xi, \eta$ and $\zeta$ correspond with the global co-ordinates, i.e. $\zeta$ coincides with the z-direction. Regarding the lower parts of Figures 12–14 it becomes obvious that—concerning the Ansatz spaces $\mathcal{S}_{\text{ts}}^{p_{\xi}, p_{\eta}, p_{\zeta}(\Omega_{\text{st}}^{h})}$ and $\mathcal{S}_{\text{ps}}^{p_{\xi}, p_{\eta}, p_{\zeta}(\Omega_{\text{st}}^{h})}$—for a wide range of p an increase of polynomial degree with respect to local co-ordinates $\xi$ and $\eta$ leads to more efficient discretizations than the ones, where the order for the local co-ordinate $\zeta$ is raised. This effect is even more evident, if we regard the Ansatz space $\mathcal{S}^{p, p, q}(\Omega_{\text{st}}^{h})$ (see lower part of Figure 14). This difference between the behaviour of the three Ansatz spaces can be explained by the anisotropic nature of the set of shape functions of $\mathcal{S}^{p, p, q}(\Omega_{\text{st}}^{h})$. While the Ansatz spaces $\mathcal{S}_{\text{ts}}^{p_{\xi}, p_{\eta}, p_{\zeta}(\Omega_{\text{st}}^{h})}$ and $\mathcal{S}_{\text{ps}}^{p_{\xi}, p_{\eta}, p_{\zeta}(\Omega_{\text{st}}^{h})}$ supply in each direction the same type of approximation, $\mathcal{S}^{p, p, q}(\Omega_{\text{st}}^{h})$ leads even for a uniform polynomial degree template p in $\zeta$-direction to a higher accuracy than in the other local directions (see Appendix A). For a wide range of p an increase of order in $\xi$ and $\eta$-direction supports for $\mathcal{S}_{\text{ts}}^{p_{\xi}, p_{\eta}, p_{\zeta}(\Omega_{\text{st}}^{h})}$ and $\mathcal{S}_{\text{ps}}^{p_{\xi}, p_{\eta}, p_{\zeta}(\Omega_{\text{st}}^{h})}$ a higher accuracy than an increase in $\zeta$-direction does. Further raise of p shows, that an increase of polynomial order in $\zeta$ becomes more important for the Ansatz spaces $\mathcal{S}_{\text{ts}}^{p_{\xi}, p_{\eta}, p_{\zeta}(\Omega_{\text{st}}^{h})}$ and $\mathcal{S}_{\text{ps}}^{p_{\xi}, p_{\eta}, p_{\zeta}(\Omega_{\text{st}}^{h})}$ because now the error in in-plane direction
<figure><img src="image_17.png" /><figcaption>Figure 14. Ansatz space $\mathcal{S}^{p, p, q}(\Omega_{\text{st}}^{h})$ with $p = 1, \dots, 6$.</figcaption></figure>

is less than that in thickness direction. Due to the built-in anisotropy, this is not the case for $\mathcal{S}^{p, p, q}(\Omega_{\text{st}}^{h})$, where the $\zeta$-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$ differently and to choose for $u_z$ a higher order approximation, e.g.:

$$ \begin{pmatrix} u_x & u_y & u_z \\ \zeta & p & p+1 \\ \eta & p & p+1 \\ \zeta & p & p+1 \end{pmatrix} $$

Furthermore the numerical investigation has shown, that it is also advantageous to choose the polynomial degree in $\xi, \eta$-direction one order higher than in $\zeta$-direction:

$$ \begin{pmatrix} u_x & u_y & u_z \\ \xi & p+1 & p+1 \\ \eta & p+1 & p+1 \\ \zeta & p & p \end{pmatrix} $$
<figure><img src="image_18.png" /><figcaption>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.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>Figure 16. Comparison of the two- and three-dimensional discretizations for a clamped plate under constant load, polynomial degrees as in Table I.</figcaption></figure>

Combining these observations, we suggest a polynomial degree template for plate-like structures being of the form

$$ \mathbf{p} = \begin{pmatrix} \xi \\ \eta \\ \xi \end{pmatrix} \begin{pmatrix} u_x & u_y & u_z \\ p_1 & p_1 & p_2 \\ p_1 & p_1 & p_2 \\ p_3 & p_3 & p_4 \end{pmatrix} $$

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 infinity. 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 difference, 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 $\mathcal{S}_{\text{ts}}^{p, \xi, p_{\eta}, p_{\zeta}}(\Omega_{\text{st}}^{h})$, $\mathcal{S}^{p, p, q}(\Omega_{\text{st}}^{h})$, $\mathcal{S}_{\text{ps}}^{p, \xi, p_{\eta}, p_{\zeta}}(\Omega_{\text{st}}^{h})$ are compared to the two-dimensional approximation of the Reissner–Mindlin problem. 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 the polynomial degree templates for the anisotropic Ansatz spaces $\mathcal{S}_{\text{ts}}^{p, \xi, p_{\eta}, p_{\zeta}}(\Omega_{\text{st}}^{h})$ and $\mathcal{S}_{\text{ps}}^{p, \xi, p_{\eta}, p_{\zeta}}(\Omega_{\text{st}}^{h})$, $\mathcal{S}^{p, p, q}(\Omega_{\text{st}}^{h})$ are listed, each row corresponding to a sequence of approximation in Figure 16. In in-plane direction ($\xi$–$\eta$-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 $\mathcal{S}_{\text{ts}}^{p, \xi, p_{\eta}, p_{\zeta}}(\Omega_{\text{st}}^{h})$ and $\mathcal{S}^{p, p, q}(\Omega_{\text{st}}^{h})$ the polynomial degree in plate thickness direction ($\zeta$) is limited for $u_x$ and $u_y$ to $p=3$ and for $u_z$ to $p=4$. Only for the Ansatz space $\mathcal{S}_{\text{ts}}^{p, \xi, p_{\eta}, p_{\zeta}}(\Omega_{\text{st}}^{h})$ the polynomial degree
<table><caption>Table I. Polynomial degree matrix for the clamped plate under uniform pressure.</caption><thead><tr><th>run</th><th>1</th><th>2</th><th>3</th><th>4</th><th>5</th><th>6</th><th>7</th><th>8</th></tr></thead><tbody><tr><td rowspan="3">S<sub>ts</sub><sup>p<sub>ζ</sub>, p<sub>η</sub>, p<sub>ζ</sub></sup>(Ω<sub>st</sub><sup>h</sup>)</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>2</td><td>2</td><td>2</td><td>3</td><td>3</td><td>3</td><td>4</td><td>4</td><td>4</td><td>5</td><td>5</td><td>5</td><td>6</td><td>6</td><td>6</td><td>7</td><td>7</td><td>7</td><td>8</td></tr><tr><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>2</td><td>2</td><td>2</td><td>3</td><td>3</td><td>3</td><td>4</td><td>4</td><td>4</td><td>5</td><td>5</td><td>5</td><td>6</td><td>6</td><td>6</td><td>7</td><td>7</td><td>7</td><td>8</td></tr><tr><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>2</td><td>2</td><td>2</td><td>3</td><td>3</td><td>3</td><td>4</td><td>4</td><td>4</td><td>5</td><td>5</td><td>5</td><td>6</td></tr><tr><td rowspan="3">S<sub>ps</sub><sup>p<sub>ζ</sub>, p<sub>η</sub>, p<sub>ζ</sub></sup>(Ω<sub>st</sub><sup>h</sup>), S<sub>p</sub><sup>p, p, q</sup>(Ω<sub>st</sub><sup>h</sup>)</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>2</td><td>2</td><td>2</td><td>3</td><td>3</td><td>3</td><td>4</td><td>4</td><td>4</td><td>5</td><td>5</td><td>5</td><td>6</td><td>6</td><td>6</td><td>7</td><td>7</td><td>7</td><td>8</td></tr><tr><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>2</td><td>2</td><td>2</td><td>3</td><td>3</td><td>3</td><td>4</td><td>4</td><td>4</td><td>5</td><td>5</td><td>5</td><td>6</td><td>6</td><td>6</td><td>7</td><td>7</td><td>7</td><td>8</td></tr><tr><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>2</td><td>2</td><td>2</td><td>3</td><td>3</td><td>3</td><td>4</td><td>3</td><td>3</td><td>4</td><td>3</td><td>3</td><td>4</td></tr></tbody></table>

<figure><img src="image_20.png" /><figcaption>Figure 17. Plate with nine columns.</figcaption></figure>

in $\zeta$-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 $\zeta$-direction is the fact, that the Ansatz space $S_{\text{ts}}^{p_{\zeta}, p_{\eta}, p_{\zeta}}(\Omega_{\text{st}}^{h})$ is not as rich as $S_{\text{ps}}^{p_{\zeta}, p_{\eta}, p_{\zeta}}(\Omega_{\text{st}}^{h})$ and $S^{p, p, q}(\Omega_{\text{st}}^{h})$ are. Therefore, $S_{\text{ts}}^{p_{\zeta}, p_{\eta}, p_{\zeta}}(\Omega_{\text{st}}^{h})$ supplies less shape functions in $\zeta$-direction and a higher polynomial degree is required to achieve approximately the same accuracy as with the Ansatz spaces $S_{\text{ps}}^{p_{\zeta}, p_{\eta}, p_{\zeta}}(\Omega_{\text{st}}^{h})$ and $S^{p, p, q}(\Omega_{\text{st}}^{h})$.

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 effort dramatically and that anisotropic $S_{\text{ts}}^{p_{\zeta}, p_{\eta}, p_{\zeta}}(\Omega_{\text{st}}^{h})$ as well as $S^{p, p, q}(\Omega_{\text{st}}^{h})$ yields a strain energy being comparable to that of the Reissner–Mindlin model with a similar number of degrees of freedom.

## 4.2. Plate with columns

As a second example we consider a plate under uniform load supported by 9 columns, where significant three-dimensional effects 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 \times 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 = 3000000$ being
<figure><img src="image_21.png" /><figcaption>Figure 18. *p*-version mesh with refinement at the boundary and at reentrant corners.</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>Figure 19. Initial mesh and adapted meshes with MITC4 elements.</figcaption></figure>

Young's modulus and $v=0.2$ denoting Poisson's ratio. In the two-dimensional case a shear correction factor $\kappa=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. Different 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 deflection *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 refined, 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 refined 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_{ps}^{p,p}(\Omega_{st}^q)$ with $p=9, 10, 11$.
<figure><img src="image_23.png" /><figcaption>Figure 20. Relative error in energy norm for the two-dimensional discretizations of a plate with columns.</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>Figure 21. Discretization of a three-dimensional plate with 194 hexahedral elements.</figcaption></figure>

<figure><img src="image_25.png" /><figcaption>Figure 22. Relative error in energy norm obtained with $(u_x, u_y) \in \mathcal{S}_{\text{ts}}^{p, p, p}(\Omega_{\text{st}}^h)$, $u_z \in \mathcal{S}_{\text{ts}}^{p+1, p+1, p+1}(\Omega_{\text{st}}^h)$.</figcaption></figure>

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 finest 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 $\mathcal{S}_{\text{ts}}^{p_2, p_0, p_2}(\Omega_{\text{st}}^h)$ is performed. The polynomial degree for the deflection $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 three-dimensional solution being estimated from a Richardson extrapolation of the strain energies
<table><caption>Table II. Computational costs of different discretizations.</caption><thead><tr><th>Run</th><th>Model</th><th>Element type</th><th>Number of elements</th><th>dof</th><th>Rel. error (%)</th><th>CPU time (s)</th></tr></thead><tbody><tr><td>1</td><td>R-M</td><td>MITC4</td><td>16978</td><td>51519</td><td>4.4</td><td>128/134</td></tr><tr><td>2</td><td>R-M</td><td>(β<sub>x</sub>, β<sub>y</sub>, w) ∈ ℓ<sub>ps</sub><sup>10,10</sup>(Ω<sub>st</sub><sup>q</sup>)</td><td>176</td><td>52950</td><td>0.9</td><td>107</td></tr><tr><td>3</td><td>3D Lamé</td><td>(u<sub>x</sub>, u<sub>y</sub>) ∈ ℓ<sub>ts</sub><sup>7,7,7</sup>(Ω<sub>st</sub><sup>h</sup>),<br>u<sub>z</sub> ∈ ℓ<sub>ts</sub><sup>8,8,8</sup>(Ω<sub>st</sub><sup>h</sup>)</td><td>194</td><td>51393</td><td>5.1</td><td>1260</td></tr></tbody></table>

obtained with a mesh consisting of 202 hexahedral elements and an Ansatz space $\mathcal{S}_{\text{ts}}^{p, p, p}(\Omega_{\text{st}}^{h})$ 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 significance 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 fine 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 three-dimensional 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 finest mesh with MITC4-elements 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 $(\beta_x, \beta_y, w) \in \mathcal{S}_{\text{ps}}^{10,10}(\Omega_{\text{st}}^{q})$ results in a computational time of only 107 s. 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) \in \mathcal{S}_{\text{ts}}^{7,7,7}(\Omega_{\text{st}}^{h})$, $u_z \in \mathcal{S}_{\text{ts}}^{8,8,8}(\Omega_{\text{st}}^{h})$ 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 significantly higher than in the two-dimensional case, it is only a small fraction of the effort which would be necessary using low-order three-dimensional elements. To efficiently integrate the (large) element stiffness matrices being the most expensive part in *p*-version computations we have further developed the vector-integration method first 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><img src="image_26.png" /><figcaption>Figure 23. Discretization of the Scordelis–Lo shell with 3 and 84 hexahedral elements.</figcaption></figure>

<figure><img src="image_27.png" /><figcaption>Figure 24. Numerical solution of the Scordelis–Lo shell with isotropic Ansatz spaces.</figcaption></figure>

<figure><img src="image_28.png" /><figcaption>Figure 25. Numerical solution of the Scordelis–Lo shell with anisotropic Ansatz spaces.</figcaption></figure>

### 4.3. Scordelis–Lo shell

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 $\rho=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 \times 10^8$ being Young’s modulus and $\nu=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 refined 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 $\mathcal{S}_{\text{ts}}^{p,p,p}(\Omega_{\text{st}}^h)$ 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$.
<figure><img src="image_29.png" /><figcaption>Figure 26. Displacement $u_y$ of point B.</figcaption></figure>

<figure><img src="image_30.png" /><figcaption>Figure 27. Cutline A–B: von Mises stress.</figcaption></figure>

In Figure 24 the relative error of computations with isotropic Ansatz spaces $\mathcal{S}_{\text{ts}}^{p_{\zeta}, p_{\eta}, p_{\zeta}}(\Omega_{\text{st}}^{\text{h}})$, $\mathcal{S}_{\text{ps}}^{p_{\zeta}, p_{\eta}, p_{\zeta}}(\Omega_{\text{st}}^{\text{h}})$ and space $\mathcal{S}^{p, p, q}(\Omega_{\text{st}}^{\text{h}})$ with $p=q=p_{\zeta}=p_{\eta}=p_{\zeta}=1, \dots, 8$ is pictured. Again, it turns out that the trunk space is the most efficient one. Further improvement of efficiency can be obtained, if anisotropic Ansatz spaces are constructed, e.g. the polynomial degree in shell thickness direction is limited to $p_{\zeta}=q=3$ (see Figure 25). If we consider e.g. the solution obtained with the trunk space $\mathcal{S}_{\text{ts}}^{p_{\zeta}, p_{\eta}, p_{\zeta}}(\Omega_{\text{st}}^{\text{h}})$ where $p_{\zeta}=p_{\eta}=8$, $p_{\zeta}=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 $\mathcal{S}_{\text{ts}}^{p_{\zeta}, p_{\eta}, p_{\zeta}}(\Omega_{\text{st}}^{\text{h}})$, $p_{\zeta}=p_{\eta}=1, \dots, 8$, $p_{\zeta}=1, \dots, 3$. In Figure 27 the von Mises stress at the cutline A–B at the middle surface of the shell is plotted. The solution, obtained with three elements and $\mathcal{S}_{\text{ts}}^{p_{\zeta}, p_{\eta}, p_{\zeta}}(\Omega_{\text{st}}^{\text{h}})$, $p_{\zeta}=p_{\eta}=8$, $p_{\zeta}=3$ with a resulting number of 594 degrees of freedom shows no deviation from the reference solution.

### 4.4. A complex shell model

Finally, we consider a more complex construction (see Figure 28), loaded by vertical weight $\rho=100$. The material behaviour is assumed to be linear elastic with $E=29000000$ being Young's modulus and $v=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 $\mathcal{S}_{\text{ts}}^{p_{\zeta}, p_{\eta}, p_{\zeta}}(\Omega_{\text{st}}^{\text{h}})$ where $p=p_{\zeta}=p_{\eta}=p_{\zeta}=1, \dots, 8$ is plotted in Figure 29.

Using a polynomial degree of $p=8$ and a corresponding number of 13408 degrees of freedom an accuracy with approximately 3.1 per cent error in energy norm is achieved. This error
<figure><img src="image_31.png" /><figcaption>Figure 28. System and mesh of a complex shell model.</figcaption></figure>

<figure><img src="image_32.png" /><figcaption>Figure 29. Relative error in energy norm of a complex shell model.</figcaption></figure>

is estimated from an extrapolation of the strain energies obtained with Ansatz space $\mathcal{S}_{\text{ts}}^{p,p,p}(\Omega_{\text{st}}^{h})$ 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 fine post-processing mesh again, being obtained by a subdivision of each *p*-element.

## 5. 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 efficiently compute the structural behaviour of plate and shell-like structures. Three different Ansatz spaces are investigated and it is shown that a flexible 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.
A. ANSATZ SPACES FOR A HEXAHEDRAL ELEMENT FORMULATION

A.1. The trunk space $\mathcal{S}_{\text{ts}}^{p_{\xi}, p_{\eta}, p_{\zeta}}(\Omega_{\text{st}}^{\text{h}})$

Eight nodal modes

$$N_{1,1,1}^{N_i}(\xi, \eta, \zeta) = \frac{1}{8}(1 + \xi\zeta)(1 + \eta\eta)(1 + \zeta\zeta), \quad i=1, \dots, 8 \qquad (\text{A1})$$

4($p_{\xi} - 1$) edge modes

4($p_{\zeta} - 1$) edge modes in $\zeta$-direction

$$\begin{align*}
N_{i,1,1}^{E_1}(\xi, \eta, \zeta) &= \frac{1}{4}(1 - \eta)(1 - \zeta)\phi_i(\xi), & i &= 2, \dots, p_{\xi} \\
N_{i,1,1}^{E_3}(\xi, \eta, \zeta) &= \frac{1}{4}(1 + \eta)(1 - \zeta)\phi_i(\xi), & i &= 2, \dots, p_{\zeta}
\end{align*} \tag{A2}$$

$$\begin{align*}
N_{i,1,1}^{E_9}(\xi, \eta, \zeta) &= \frac{1}{4}(1 - \eta)(1 + \zeta)\phi_i(\xi), & i &= 2, \dots, p_{\zeta} \\
N_{i,1,1}^{E_{11}}(\xi, \eta, \zeta) &= \frac{1}{4}(1 + \eta)(1 + \zeta)\phi_i(\xi), & i &= 2, \dots, p_{\xi}
\end{align*}$$

4($p_{\eta} - 1$) edge modes in $\eta$-direction

$$\begin{align*}
N_{1,j,1}^{E_2}(\xi, \eta, \zeta) &= \frac{1}{4}(1 + \xi)(1 - \zeta)\phi_j(\eta), & j &= 2, \dots, p_{\eta} \\
N_{1,j,1}^{E_4}(\xi, \eta, \zeta) &= \frac{1}{4}(1 - \xi)(1 - \zeta)\phi_j(\eta), & j &= 2, \dots, p_{\zeta}
\end{align*} \tag{A3}$$

$$\begin{align*}
N_{1,j,1}^{E_{10}}(\xi, \eta, \zeta) &= \frac{1}{4}(1 + \xi)(1 + \zeta)\phi_j(\eta), & j &= 2, \dots, p_{\eta} \\
N_{1,j,1}^{E_{12}}(\xi, \eta, \zeta) &= \frac{1}{4}(1 - \xi)(1 + \zeta)\phi_j(\eta), & j &= 2, \dots, p_{\zeta}
\end{align*}$$

4($p_{\zeta} - 1$) edge modes in $\zeta$-direction

$$\begin{align*}
N_{1,1,k}^{E_5}(\xi, \eta, \zeta) &= \frac{1}{4}(1 - \xi)(1 - \eta)\phi_k(\zeta), & k &= 2, \dots, p_{\zeta} \\
N_{1,1,k}^{E_6}(\xi, \eta, \zeta) &= \frac{1}{4}(1 + \xi)(1 - \eta)\phi_k(\zeta), & k &= 2, \dots, p_{\zeta}
\end{align*} \tag{A4}$$

$$\begin{align*}
N_{1,1,k}^{E_7}(\xi, \eta, \zeta) &= \frac{1}{4}(1 + \xi)(1 + \eta)\phi_k(\zeta), & k &= 2, \dots, p_{\zeta} \\
N_{1,1,k}^{E_8}(\xi, \eta, \zeta) &= \frac{1}{4}(1 - \xi)(1 + \eta)\phi_k(\zeta), & k &= 2, \dots, p_{\zeta}
\end{align*}$$
Face modes

$$
N_{i,j,1}^{F_1}(\xi, \eta, \zeta) = \frac{1}{2}(1-\zeta)\phi_i(\xi)\phi_j(\eta), \quad i=2, \dots, p_{\xi}-2, \quad j=2, \dots, p_{\eta}-2, \quad i+j=4, \dots, \max\{p_{\xi}, p_{\eta}\}
$$

$$
N_{i,1,k}^{F_2}(\xi, \eta, \zeta) = \frac{1}{2}(1-\eta)\phi_i(\xi)\phi_k(\zeta), \quad i=2, \dots, p_{\xi}-2, \quad k=2, \dots, p_{\zeta}-2, \quad i+k=4, \dots, \max\{p_{\xi}, p_{\zeta}\}
$$

$$
N_{1,j,k}^{F_3}(\xi, \eta, \zeta) = \frac{1}{2}(1+\zeta)\phi_j(\eta)\phi_k(\zeta), \quad j=2, \dots, p_{\eta}-2, \quad k=2, \dots, p_{\zeta}-2, \quad j+k=4, \dots, \max\{p_{\eta}, p_{\zeta}\}
$$

$$
N_{i,1,k}^{F_4}(\xi, \eta, \zeta) = \frac{1}{2}(1+\eta)\phi_i(\xi)\phi_k(\zeta), \quad i=2, \dots, p_{\xi}-2, \quad k=2, \dots, p_{\zeta}-2, \quad i+k=4, \dots, \max\{p_{\xi}, p_{\zeta}\}
$$

$$
N_{1,j,k}^{F_5}(\xi, \eta, \zeta) = \frac{1}{2}(1-\xi)\phi_j(\eta)\phi_k(\zeta), \quad j=2, \dots, p_{\eta}-2, \quad k=2, \dots, p_{\zeta}-2, \quad j+k=4, \dots, \max\{p_{\eta}, p_{\zeta}\}
$$

$$
N_{i,j,1}^{F_6}(\xi, \eta, \zeta) = \frac{1}{2}(1+\zeta)\phi_i(\xi)\phi_j(\eta), \quad i=2, \dots, p_{\xi}-2, \quad j=2, \dots, p_{\eta}-2, \quad i+j=4, \dots, \max\{p_{\xi}, p_{\eta}\}
$$

(A5)

Internal modes

$$
N_{i,j,k}^{\text{int}}(\xi, \eta, \zeta) = \phi_i(\xi)\phi_j(\eta)\phi_k(\zeta) \quad i=2, \dots, p_{\xi}-4, \quad j=2, \dots, p_{\eta}-4, \quad k=2, \dots, p_{\zeta}-4,
$$

$$
i + j + k = 6, \dots, \max\{p_{\xi}, p_{\eta}, p_{\zeta}\} \tag{A6}
$$

A.2. The tensor product space $\mathcal{S}_{\text{ps}}^{p_{\xi}, p_{\eta}, p_{\zeta}}(\Omega_{\text{st}}^{h})$

Nodal and edge modes as in Equations (19)–(22)

$$
2[(p_{\xi}-1)(p_{\eta}-1) + (p_{\xi}-1)(p_{\zeta}-1) + (p_{\eta}-1)(p_{\zeta}-1)] \text{ face modes}
$$

$$
N_{i,j,1}^{F_1}(\xi, \eta, \zeta) = \frac{1}{2}(1-\zeta)\phi_i(\xi)\phi_j(\eta), \quad i=2, \dots, p_{\xi}, \quad j=2, \dots, p_{\eta}
$$

$$
N_{i,1,k}^{F_2}(\xi, \eta, \zeta) = \frac{1}{2}(1-\eta)\phi_i(\xi)\phi_k(\zeta), \quad i=2, \dots, p_{\xi}, \quad k=2, \dots, p_{\zeta}
$$

$$
N_{1,j,k}^{F_3}(\xi, \eta, \zeta) = \frac{1}{2}(1+\xi)\phi_j(\eta)\phi_k(\zeta), \quad j=2, \dots, p_{\eta}, \quad k=2, \dots, p_{\zeta}
$$

(A7)

$$
N_{i,1,k}^{F_4}(\xi, \eta, \zeta) = \frac{1}{2}(1+\eta)\phi_i(\xi)\phi_k(\zeta), \quad i=2, \dots, p_{\xi}, \quad k=2, \dots, p_{\zeta}
$$

$$
N_{1,j,k}^{F_5}(\xi, \eta, \zeta) = \frac{1}{2}(1-\xi)\phi_j(\eta)\phi_k(\zeta), \quad j=2, \dots, p_{\eta}, \quad k=2, \dots, p_{\zeta}
$$

$$
N_{i,j,1}^{F_6}(\xi, \eta, \zeta) = \frac{1}{2}(1+\zeta)\phi_i(\xi)\phi_j(\eta), \quad i=2, \dots, p_{\xi}, \quad j=2, \dots, p_{\eta}
$$
$(p_{\xi}-1)(p_{\eta}-1)(p_{\zeta}-1) \text{ internal modes}$

$$
N_{i,j,k}^{\text{int}}(\xi, \eta, \zeta) = \phi_i(\xi)\phi_j(\eta)\phi_k(\zeta), \quad i=2, \dots, p_{\xi}, \quad j=2, \dots, p_{\eta}, \quad k=2, \dots, p_{\zeta} \quad (\text{A8})
$$

A.3. The anisotropic tensor product space $\mathcal{S}^{p, p, q}(\Omega_{st}^{h})$

Nodal modes as in (19)

$8(p-1)+4(q-1)$ edge modes

8( p − 1) edge modes in ξ − η plane

$$
\begin{align*}
N_{i,1,1}^{E_1}(\xi, \eta, \zeta) &= \frac{1}{4}(1-\eta)(1-\zeta)\phi_i(\xi) & i=2, \dots, p \\
N_{1,j,1}^{E_2}(\xi, \eta, \zeta) &= \frac{1}{4}(1+\xi)(1-\zeta)\phi_j(\eta) & j=2, \dots, p \\
N_{i,1,1}^{E_3}(\xi, \eta, \zeta) &= \frac{1}{4}(1+\eta)(1-\zeta)\phi_i(\xi) & i=2, \dots, p \\
N_{1,j,1}^{E_4}(\xi, \eta, \zeta) &= \frac{1}{4}(1-\xi)(1-\zeta)\phi_j(\eta) & j=2, \dots, p \\
N_{i,1,1}^{E_9}(\xi, \eta, \zeta) &= \frac{1}{4}(1-\eta)(1+\zeta)\phi_i(\xi) & i=2, \dots, p \\
N_{1,j,1}^{E_{10}}(\xi, \eta, \zeta) &= \frac{1}{4}(1+\xi)(1+\zeta)\phi_j(\eta) & j=2, \dots, p \\
N_{i,1,1}^{E_{11}}(\xi, \eta, \zeta) &= \frac{1}{4}(1+\eta)(1+\zeta)\phi_i(\xi) & i=2, \dots, p \\
N_{1,j,1}^{E_{12}}(\xi, \eta, \zeta) &= \frac{1}{4}(1-\xi)(1+\zeta)\phi_j(\eta) & j=2, \dots, p
\tag{A9}
\end{align*}
$$

4( q − 1) edge modes in ζ -direction

$$
\begin{align*}
N_{1,1,k}^{E_5}(\xi, \eta, \zeta) &= \frac{1}{4}(1-\xi)(1-\eta)\phi_k(\zeta) & k=2, \dots, q \\
N_{1,1,k}^{E_6}(\xi, \eta, \zeta) &= \frac{1}{4}(1+\xi)(1-\eta)\phi_k(\zeta) & k=2, \dots, q \\
N_{1,1,k}^{E_7}(\xi, \eta, \zeta) &= \frac{1}{4}(1+\xi)(1+\eta)\phi_k(\zeta) & k=2, \dots, q \\
N_{1,1,k}^{E_8}(\xi, \eta, \zeta) &= \frac{1}{4}(1-\xi)(1+\eta)\phi_k(\zeta) & k=2, \dots, q
\tag{A10}
\end{align*}
$$

$4(p-1)(q-1)+(p-2)(p-3)$ face modes
face modes in $\xi-\eta$ plane

$$
\begin{align}
N_{i,j,1}^{\mathrm{F}_1}(\xi, \eta, \zeta) &= \frac{1}{2}(1-\zeta)\phi_i(\xi)\phi_j(\eta), \quad i,j=2, \dots, p-2, \quad i+j=4, \dots, p \tag{A11} \\
N_{i,j,1}^{\mathrm{F}_6}(\xi, \eta, \zeta) &= \frac{1}{2}(1+\zeta)\phi_i(\xi)\phi_j(\eta), \quad i,j=2, \dots, p-2, \quad i+j=4, \dots, p \nonumber
\end{align}
$$
<figure><img src="image_33.png" /><figcaption>Plate 1. Deformation (scaling factor = 100) of a three-dimensional plate with columns.</figcaption></figure>

<figure><img src="image_34.png" /><figcaption>Plate 2. Deformation (scaling factor = 500) and von Mises stress of a complex shell model with a zoomed detail.</figcaption></figure>
face modes in $\xi-\zeta$ and $\eta-\zeta$ plane

$$
\begin{align*}
N_{i,1,k}^{F_2}(\xi, \eta, \zeta) &= \frac{1}{2}(1-\eta)\phi_i(\xi)\phi_k(\zeta), & i=2,\dots,p, \quad k=2,\dots,q \\
N_{1,j,k}^{F_3}(\xi, \eta, \zeta) &= \frac{1}{2}(1+\xi)\phi_j(\eta)\phi_k(\zeta), & j=2,\dots,p, \quad k=2,\dots,q \\
N_{i,1,k}^{F_4}(\xi, \eta, \zeta) &= \frac{1}{2}(1+\eta)\phi_i(\xi)\phi_k(\zeta), & i=2,\dots,p, \quad k=2,\dots,q \\
N_{1,j,k}^{F_5}(\xi, \eta, \zeta) &= \frac{1}{2}(1-\xi)\phi_j(\eta)\phi_k(\zeta), & j=2,\dots,p, \quad k=2,\dots,q
\end{align*}
\tag{A12}
$$

$(p-3)(p-2)(q-1)/2$ internal modes

$$
N_{i,j,k}^{\text{int}}(\xi, \eta, \zeta) = \phi_i(\xi)\phi_j(\eta)\phi_k(\zeta), \quad i,j=2,\dots,p-2, \quad i+j=4,\dots,p, \quad k=2,\dots,q \quad \text{(A13)}
$$

B. THE BLENDING FUNCTION METHOD FOR HEXAHEDRAL ELEMENTS

B.1. Edge blending (from Reference [15])

$$
\begin{align*}
\text{Edge 1: } \mathbf{e}_1(\xi, \eta, \zeta) &:= \left( \mathbf{E}_1(\xi) - \frac{(1-\xi)\mathbf{X}_1 + (1+\xi)\mathbf{X}_2}{2} \right) \left( \frac{1-\eta}{2} \right) \left( \frac{1-\zeta}{2} \right) \\
\text{Edge 2: } \mathbf{e}_2(\xi, \eta, \zeta) &:= \left( \mathbf{E}_2(\eta) - \frac{(1-\eta)\mathbf{X}_2 + (1+\eta)\mathbf{X}_3}{2} \right) \left( \frac{1+\xi}{2} \right) \left( \frac{1-\zeta}{2} \right) \\
\text{Edge 3: } \mathbf{e}_3(\xi, \eta, \zeta) &:= \left( \mathbf{E}_3(\xi) - \frac{(1-\xi)\mathbf{X}_4 + (1+\xi)\mathbf{X}_5}{2} \right) \left( \frac{1+\eta}{2} \right) \left( \frac{1-\zeta}{2} \right) \\
\text{Edge 4: } \mathbf{e}_4(\xi, \eta, \zeta) &:= \left( \mathbf{E}_4(\eta) - \frac{(1-\eta)\mathbf{X}_1 + (1+\eta)\mathbf{X}_4}{2} \right) \left( \frac{1-\zeta}{2} \right) \left( \frac{1-\zeta}{2} \right) \\
\text{Edge 5: } \mathbf{e}_5(\xi, \eta, \zeta) &:= \left( \mathbf{E}_5(\zeta) - \frac{(1-\zeta)\mathbf{X}_1 + (1+\zeta)\mathbf{X}_5}{2} \right) \left( \frac{1-\xi}{2} \right) \left( \frac{1-\eta}{2} \right) \tag{B1} \\
\text{Edge 6: } \mathbf{e}_6(\xi, \eta, \zeta) &:= \left( \mathbf{E}_6(\zeta) - \frac{(1-\zeta)\mathbf{X}_2 + (1+\zeta)\mathbf{X}_6}{2} \right) \left( \frac{1+\xi}{2} \right) \left( \frac{1-\eta}{2} \right) \\
\text{Edge 7: } \mathbf{e}_7(\xi, \eta, \zeta) &:= \left( \mathbf{E}_7(\zeta) - \frac{(1-\zeta)\mathbf{X}_3 + (1+\zeta)\mathbf{X}_7}{2} \right) \left( \frac{1+\xi}{2} \right) \left( \frac{1+\eta}{2} \right) \\
\text{Edge 8: } \mathbf{e}_8(\xi, \eta, \zeta) &:= \left( \mathbf{E}_8(\zeta) - \frac{(1-\zeta)\mathbf{X}_4 + (1+\zeta)\mathbf{X}_8}{2} \right) \left( \frac{1-\xi}{2} \right) \left( \frac{1+\eta}{2} \right) \\
\text{Edge 9: } \mathbf{e}_9(\xi, \eta, \zeta) &:= \left( \mathbf{E}_9(\xi) - \frac{(1-\xi)\mathbf{X}_5 + (1+\xi)\mathbf{X}_6}{2} \right) \left( \frac{1-\eta}{2} \right) \left( \frac{1+\zeta}{2} \right)
\end{align*}
$$
$$ \text{Edge 10: } \mathbf{e}_{10}(\xi, \eta, \zeta) := \left( \mathbf{E}_{10}(\eta) - \frac{(1-\eta)\mathbf{X}_6 + (1+\eta)\mathbf{X}_7}{2} \right) \left( \frac{1+\xi}{2} \right) \left( \frac{1-\xi}{2} \right) $$

$$ \text{Edge 11: } \mathbf{e}_{11}(\xi, \eta, \zeta) := \left( \mathbf{E}_{11}(\xi) - \frac{(1-\xi)\mathbf{X}_8 + (1+\xi)\mathbf{X}_7}{2} \right) \left( \frac{1+\eta}{2} \right) \left( \frac{1-\zeta}{2} \right) $$

$$ \text{Edge 12: } \mathbf{e}_{12}(\xi, \eta, \zeta) := \left( \mathbf{E}_{12}(\eta) - \frac{(1-\eta)\mathbf{X}_5 + (1+\eta)\mathbf{X}_8}{2} \right) \left( \frac{1-\zeta}{2} \right) \left( \frac{1+\zeta}{2} \right) $$

## B.2. Face blending (from Reference [15])

$$ \text{Face 1: } \mathbf{f}_1(\xi, \eta, \zeta) := (\mathbf{F}_1(\xi, \eta) - \frac{1}{4}((1-\xi)(1-\eta)\mathbf{X}_1 + (1+\xi)(1-\eta)\mathbf{X}_2 + (1+\xi)(1+\eta)\mathbf{X}_3 + (1-\xi)(1+\eta)\mathbf{X}_4)) \left(\frac{1-\zeta}{2}\right) $$

$$ \text{Face 2: } \mathbf{f}_2(\xi, \eta, \zeta) := (\mathbf{F}_2(\xi, \zeta) - \frac{1}{4}((1-\xi)(1-\zeta)\mathbf{X}_1 + (1+\xi)(1-\zeta)\mathbf{X}_2 + (1+\xi)(1+\zeta)\mathbf{X}_3 + (1-\xi)(1+\zeta)\mathbf{X}_5)) \left(\frac{1-\eta}{2}\right) $$

$$ \text{Face 3: } \mathbf{f}_3(\xi, \eta, \zeta) := (\mathbf{F}_3(\eta, \zeta) - \frac{1}{4}((1-\eta)(1-\zeta)\mathbf{X}_2 + (1+\eta)(1-\zeta)\mathbf{X}_3 + (1+\eta)(1+\zeta)\mathbf{X}_7 + (1-\eta)(1+\zeta)\mathbf{X}_6)) \left(\frac{1+\zeta}{2}\right) \quad (\text{B2}) $$

$$ \text{Face 4: } \mathbf{f}_4(\xi, \eta, \zeta) := (\mathbf{F}_4(\xi, \zeta) - \frac{1}{4}((1-\xi)(1-\zeta)\mathbf{X}_4 + (1+\xi)(1-\zeta)\mathbf{X}_3 + (1+\xi)(1+\zeta)\mathbf{X}_7 + (1-\xi)(1+\zeta)\mathbf{X}_8)) \left(\frac{1+\eta}{2}\right) $$

$$ \text{Face 5: } \mathbf{f}_5(\xi, \eta, \zeta) := (\mathbf{F}_5(\eta, \zeta) - \frac{1}{4}((1-\eta)(1-\zeta)\mathbf{X}_1 + (1+\eta)(1-\zeta)\mathbf{X}_4 + (1+\eta)(1+\zeta)\mathbf{X}_8 + (1-\eta)(1+\zeta)\mathbf{X}_5)) \left(\frac{1-\zeta}{2}\right) $$

$$ \text{Face 6: } \mathbf{f}_6(\xi, \eta, \zeta) := (\mathbf{F}_6(\xi, \eta) - \frac{1}{4}((1-\xi)(1-\eta)\mathbf{X}_5 + (1+\xi)(1-\eta)\mathbf{X}_6 + (1+\xi)(1+\eta)\mathbf{X}_7 + (1-\xi)(1+\eta)\mathbf{X}_8)) \left(\frac{1+\zeta}{2}\right) $$

## REFERENCES

1. Babuška I, Szabó B. On the rates of convergence of the finite element method. *International Journal for Numerical Methods in Engineering* 1982; **18**:323–341.

2. Holzer S, Rank E, Werner H. An implementation of the hp-version of the finite element method for Reissner-Mindlin plate problems. *International Journal for Numerical Methods in Engineering* 1990; **30**:459–471.

3. Rank E, Krause R, Preusch K. On the accuracy of p-version elements for the Reissner-Mindlin plate problem. *International Journal for Numerical Methods in Engineering* 1998; **43**:51–67.

4. Rank E, Düster A, Krafczyk M, Rücker M. Some aspects of coupling structural models and p-version finite element models. *Proceedings of the 4th World Congress on Computational Mechanics*, Buenos Aires, Argentina, 1998.
5. Rank E, Rücker M, Düster A, Bröker H. The efficiency of the *p*-version finite element method. *Proceedings of ECCM '99, European Conference on Computational Mechanics*, München, Germany, CD-ROM, 1999.

6. Schwab C. *p- and hp-finite element methods, theory and applications in solid and fluid mechanics*. Oxford University Press: Oxford, 1998.

7. Szabó B, Babuška I. *Finite element analysis*. Wiley: New York, 1991.

8. Szabó B. Mesh design for the *p*-version of the finite element method. *Computer Methods in Applied Mechanics and Engineering* 1986; **55**:181–197.

9. Suri M. Analytical and computational assessment of locking in the hp finite element method. *Computer Methods in Applied Mechanics and Engineering* 1996; **133**:347–371.

10. Düster A, Rank E. The *p*-version of the finite element method compared to an adaptive *h*-version for the deformation theory of plasticity. *Computer Methods in Applied Mechanics and Engineering* 2001; **190**: 1925–1935.

11. Holzer S, Yosibash Z. The *p*-version of the finite element method in incremental elasto-plastic analysis. *International Journal for Numerical Methods in Engineering* 1996; **39**:1859–1878.

12. Szabó B, Actis R, Holzer S. Solution of elastic-plastic stress analysis problems by the *p*-version of the finite element method. In *Modeling, Mesh Generation, and Adaptive Numerical Methods for Partial Differential Equations*, IMA Volumes in Mathematics and its Applications, Babuška I, Flaherty J et al. (eds). vol. 75. Springer: New York, 1995; 395–416.

13. Actis RL, Szabó B, Schwab C. Hierarchic models for laminated plates and shells. *Computer Methods in Applied Mechanics and Engineering* 1999; **172**:79–107.

14. Gordon WJ, Hall CA. Construction of curvilinear co-ordinate systems and applications to mesh generation. *International Journal for Numerical Methods in Engineering* 1973; **7**:461–477.

15. Királyfalvi G, Szabó B. Quasi-regional mapping for the *p*-version of the finite element method. *Finite elements in analysis and design* 1997; **27**:85–97.

16. Dey S. Geometry-based three-dimensional hp-finite element modelling and computations. *Ph.D. Thesis*, Rensselaer Polytechnic Institute Troy, New York, 1997.

17. Dey S, Shephard MS, Flaherty JE. Geometry representation issues associated with *p*-version finite element computations. *Computer Methods in Applied Mechanics and Engineering* 1997; **150**:39–55.

18. Zienkiewicz OC, Craig A. Adaptive refinement, error estimates, multigrid solution, and hierarchic finite element method concepts. In *Accuracy Estimates and Adaptive Refinements in Finite Element Computations*, Babuška I, Zienkiewicz OC, Gago J, Olivera ER de A (eds). Wiley: New York, 1986.

19. Ainsworth M. A preconditioner based on domain decomposition for hp-fe approximation on quasi-uniform meshes. *SIAM Journal on Numerical Analysis* 1996; **33**(4):1358–1376.

20. Mandel J. Iterative Solvers by substructuring for the *p*-version finite element method. *Computer Methods in Applied Mechanics and Engineering* 1990; **80**:117–128.

21. Papadrakakis M, Babilis GB. Solution techniques for the *p*-version of the adaptive fe method. *International Journal for Numerical Methods in Engineering* 1994; **37**:1413–1431.

22. Bröker H. Integration von geometrischer Modellierung und Berechnung nach der *p*-Version der FEM. *Ph.D. Thesis*, Technische Universität München, 2001.

23. Stein E, Barthold F-J, Ohnimus S, Schmidt M. Adaptive Finite Elements in Elastoplasticity with mechanical Error Indicators and Neumann-type Estimators. *Proceedings of the Workshop on New Advances in Adaptive Computational Mechanics*, Cachan, September 1997.

24. Bathe K-J, Dvorkin EN. Short communication: a four-node plate bending element based on Mindlin/Reissner plate theory and a mixed interpolation. *International Journal for Numerical Methods in Engineering* 1985; **21**:367–383.

25. Düster A. Die hp-d Methode für Reissner–Mindlin Plattenprobleme. *Forum Bauinformatik, Dresden '97, Fortschritt–Berichte*, Reihe 4, No. 140. VDI-Verlag: Düsseldorf, 1997.

26. Hinnant HE. A fast method of numerical quadrature for *p*-version finite element matrices. *International Journal for Numerical Methods in Engineering* 1994; **37**:3723–3750.

27. Nübel V, Düster A, Rank E. Die Methode Vektorintegration bei der *p*-Version der Finite-Elemente-Methode. *Technical Report*, Lehrstuhl für Bauinformatik, Technische Universität München, 2000.

28. Nübel V, Düster A, Rank E. Adaptive vector integration as an efficient quadrature scheme for *p*-version finite element matrices. *European Conference on Computational Mechanics*, Cracow, Poland, 2001, submitted for publication.

29. Mandel J. Iterative methods for *p*-version finite element method: preconditioning thin solids. *Computer Methods in Applied Mechanics and Engineering* 1996; **133**:247–257.

30. Scordelis AC, Lo KS. Computer analysis of cylindrical shells. *Journal of American Concrete Institute* 1969; **61**:539–561.

31. Corney J. *3D Modeling with the ACIS Kernel and Toolkit*. Wiley: Chister, 1997.