Convergence in the incompressible limit of new discontinuous Galerkin methods with general quadrilateral and hexahedral elements

Standard low-order finite elements, which perform well for problems involving compressible elastic materials, are known to under-perform when nearly incompressible materials are involved, commonly exhibiting the locking phenomenon. Interior penalty (IP) discontinuous Galerkin methods have been shown to circumvent locking when simplicial elements are used. The same IP methods, however, result in locking on meshes of quadrilaterals. The authors have shown in earlier work that under-integration of specified terms in the IP formulation eliminates the locking problem for rectangular elements. Here it is demonstrated through an extensive numerical investigation that the effect of using under-integration carries over successfully to meshes of more general quadrilateral elements, as would likely be used in practical applications, and results in accurate displacement approximations. Uniform convergence with respect to the compressibility parameter is shown numerically. Additionally, a stress approximation obtained here by postprocessing shows good convergence in the incompressible limit.

certainly for linear problems, and for an increasingly wide range of nonlinear problems. Despite this, significant challenges remain.
In the context of solid mechanics, and in particular in problems for elastic materials, the standard Galerkin (SG) finite element method, while performing very well for compressible materials, may exhibit the phenomenon known as "volumetric locking" for materials that are nearly incompressible, if low-order (linear or bilinear, or trilinear in three dimensions) elements are used. Manifesting particularly in the case of bending-dominated problems, this pathological behaviour results from the too-severe constraint placed on the solution by the incompressibility condition. The adverse effect of the degree of compressibility on the performance of the SG method may nevertheless also be seen in poor displacement approximations that are not specifically of the locking type.
The problem may be circumvented by the use of high-order elements. Low-order elements remain an attractive option, though, and for this reason various extensions to the SG method have been studied, and shown to be effective in remedying locking when low-order approximations for the displacement are used.
One class of extensions is mixed methods (see, for example, [6]). Related to the pressure-displacement mixed method is the method of selective reduced integration (SRI), also effective in producing lockingfree results. Finally, discontinuous Galerkin (DG) methods, specifically the range of interior penalty (IP) DG methods, have been used effectively with low-order elements, within a limited scope.
Various DG mixed methods have been proven to be robust for near-incompressible isotropic elasticity, allowing overall for meshes of quadrilateral, hexahedral, as well as simplicial elements (see [15,9,17,25]).
In contrast, DG primal formulations have been established as having optimal performance independent of material parameters for meshes of triangular elements ( [23,24,15,16,7]), but not for meshes of quadrilateral elements.
In a numerical investigation, Liu et al. [18] consider the matter of locking in the context of low-order hexahedral elements. This work, on a specific benchmark problem, shows the under-performance of the SG method and the superiority of three well-known primal IP methods, NIPG, SIPG and IIPG (Nonsymmetric, Symmetric and Incomplete Interior Penalty Galerkin methods, developed in [21,20,22,11]), as well as of the method of Oden, Babuska and Baumann [19] (known as OBB), a penalty-free version of NIPG. However, the conditions of the benchmark problem are not those leading to the severest form of locking, and no accompanying analysis or convergence data is included. Therefore, while the results are positive, questions remain about the scope of the superior performance of the IP methods with non-simplicial elements.
In [13], the authors showed using several numerical examples in two and three dimensions that the three IP methods do in fact produce poor approximations, and notably locking-type behaviour, when quadrilateral/hexahedral elements are used. A new method, in which selected edge terms of the IP formulation are under-integrated, circumvents the problem, as shown through an analytical proof that the new method is locking-free for rectangular elements, and through numerical examples that demonstrate the optimal performance of this formulation.
While the technique of under-integration (or SRI) has long been used in other contexts to eliminate volume locking, it has until recently typically been used on integrals on element domains, while in the formulation of [13] it is used in integrals on element edges only. Hansbo and Larson [16] use underintegration similarly on an edge-based stabilization term in their formulation for triangular elements, relaxing what would be otherwise be a severely constraining term. In other applications, edge underintegration is used in [14] with simplicial elements to eliminate extensional locking within transversely isotropic elasticity, and in [3] this technique is used to circumvent shear locking in beams.
The analysis in [13] of the new IP formulation accommodates both essential and natural boundary conditions, and the numerical examples incorporate both, as would be expected in a realistic, practical model. However, both the theoretical and numerical analyses presented are concerned with rectangular elements and simple domains. A single example in that paper shows the method locking-free with quadrilaterals as well, but is not conclusive regarding the general case. In a more recent computational paper, Bayat et al. [3] have similarly shown several variants of IIPG to be volumetrically locking-free when selected edge terms are under-integrated. However, these authors consider limited test cases only and their results are therefore inconclusive regarding the general case.
For the method of [13] to be of practical value in solving a broad range of boundary value problems, it is necessary to establish its robustness when general quadrilaterals or hexahedrals are employed, allowing for its use on a variety of domain shapes and meshes (particularly unstructured).
This paper seeks to address this issue computationally through an extensive numerical investigation using non-rectangular elements with a variety of model problems. We systematically compare the performance of the new method, particularly in the incompressible limit, as meshes with decreasing element shape regularity are used. We demonstrate that the formulation of [13] is effective in alleviating locking and generating accurate approximations in the general case.
In practice the stress field generated from the displacement approximation by postprocessing is often of interest. As a second component of this work we therefore study the accuracy of the stress field approximation obtained from the new IP methods, considering both error convergence rates and approximation quality at individual refinement levels.
Following this introduction, the boundary value problem and DG framework and formulation are presented in §2. Section 3 gives the numerical results for both the displacement and the stress approximations, before the conclusion in §4.

The boundary value problem of linear elasticity
Let a homogeneous, isotropic, linear elastic body occupy the bounded domain Ω ⊂ R d (d = 2, 3), with the Lipschitz boundary ∂Ω consisting of a Dirichlet portion, Γ D , of positive measure, and a Neumann portion, Γ N , such that Γ D ∩ Γ N = ∅ and Γ D ∪ Γ N = ∂Ω, and with outward unit normal n.
The resultant displacement is u, and the strain ε is expressed as a tensor defined in index notation as The stress σ (u) is related to the strain via the constitutive law where λ and the shear modulus µ are known as the Lamé parameters, and 1 is the second-order identity The governing equation of the system is the equilibrium equation and the boundary conditions are The Lamé parameters λ and µ are assumed to be positive, and can be expressed in terms of the Young's modulus, E, and Poisson's ratio, ν, by As ν → 1 2 , which corresponds to the incompressible limit, so λ → ∞.

The discontinuous Galerkin framework
The outward unit normal of Ω e is denoted by n e .
In the following, all definitions and notation are given for d = 2, but are equally applicable to d = 3 if "edge" is replaced with "face" in each instance.
Each element has a boundary ∂Ω e , consisting of edges E. Define h E : = diam (E).
The union of all edges lying in the interior of the domain, rather than on the boundary, will be denoted by Γ int . Define Γ iD : = Γ int ∪ Γ D . By abuse of notation, any symbol denoting a union of edges will also denote the corresponding set of edges.
Use will also be made of the discrete Sobolev space Some kind of weak continuity is nevertheless required between neighbouring elements. For a vector v and a tensor τ , with components in H 1 (Ω e ) and H 1 (Ω f ) for adjacent elements Ω e and Ω f with common edge Γ ef , the jumps are defined by v : = v e ⊗ n e + v f ⊗ n f , τ : = τ e n e + τ f n f and [[v]] : = v e · n e + v f · n f and the averages by On edges E such that E ∩ ∂Ω = ∅, the jumps and averages are defined by v = v ⊗ n, τ = τ n, With Q 1 (Ω) the space of polynomials on Ω with maximum degree one in each variable, define the DG solution space

Modified interior penalty (IP) formulation
In the general formulation presented in [13], which defines three IP methods, selected edge terms are under-integrated (indicated below by a summation over the Gauss points, p i , for i = 1, ..., ngp, with w i referring to the corresponding weighting value). With the context here being the use of bilinear or trilinear elements, the appropriate under-integration employs a single Gauss point (that is, ngp = 1).
With non-negative parameters k µ and k λ , and a switch θ to distinguish between methods, the formulation is defined by the bilinear form and linear functional The standard IP formulation is recovered by integrating fully all terms (i.e. reverting to ngp = 2), in which case θ = 1 gives the NIPG method, θ = −1 gives SIPG, and θ = 0 gives IIPG. (The original methods do not in all cases contain the stabilization terms included here.) As in the original IP methods, both Dirichlet and Neumann boundary conditions, (2b) and (2c), are imposed weakly through this formulation.
The formulation of (6) is shown in [13] to be stable and optimally convergent, and specifically uniformly convergent in the incompressible limit, provided that the domain, applied forces and boundary conditions satisfy the necessary smoothness requirements (as detailed in [13]). In order to ascertain the effects of deviating from the use of rectangular elements with the IP methods, in the first, second and fourth examples, comparative performance of the methods is investigated on meshes with systematically increasing degrees of mesh distortion. At each refinement level, an initial mesh is generated by isotropically refining from a single element. The resulting elements are squares in two dimensions or cubes in three dimensions. Meshes are then modified according to an algorithm (applied within the deal.ii library [2]) that distorts each element up to a prescribed distortion factor (df ), higher values of df producing meshes with lower element shape regularity. Meshes of df = 0.0 (initial, isotropically refined meshes), df = 0.1 and df = 0.3 are considered here. (See Figure 1 for an example at a given refinement level.)

Computational examples
In the third example, the L-shaped domain, an unstructured, graded initial mesh is used, with a wide range of element regularity included. This mesh is then refined to investigate convergence behaviour. It should be noted that in general the shape regularity of the elements will be increased by the refinement process in this example.
Unless otherwise indicated, the IP stabilization parameters are set at k µ = k λ = 10 for the IIPG and SIPG methods, and k µ = 10, k λ = 0 for the NIPG method.
Results are given for both the original (with full edge integration) and the new IP methods. It may be assumed henceforth that references to the IP method are to the new method of [13], unless explicit mention is made of the original IP methods.
Results obtained using the SG method with Q 1 and Q 2 elements, and Q 1 with SRI, have been included for perspective on the effects of deviating from rectangular elements in both the compressible and near-incompressible cases.
Convergence plots of the H 1 error for displacement approximations and the L 2 error for post-processed stresses are displayed, where the mesh measure h is the average element diagonal. Contour plots for the n no. els no. dofs SG Q 1 SG Q 2  IP  1  4  18  50  32  2  16  50  162  128  3  64  162  578  512  4  256  578  2178  2048  5  1024  2178  8450  8192  6  4096  8450  33282  32768  7 16384  33282  132098  131072  8  Where the exact solution is shown in a contour plot, nodal values are calculated directly from the analytical solution at nodal points, for both displacement and stress fields, except in the case of the L-shaped domain. There, for the stress field, the exact values are calculated from the analytical solution at quadrature points, and a projection, like that for the approximate stresses, is done to obtain nodal values. The reason and implications will be discussed where the results of that example are given.

Cantilever beam
A square beam in two dimensions, with E = 1500000, is subjected to a linearly varying force on the free end, with the maximum value f = 3000, as illustrated in Figure 2. The analytical solution is given in [12].
Eight levels of mesh refinement are used, and Table 1 details the number of elements and degrees of freedom (dofs) for each method at each level.
Displacement results (scaled) of the original and new IP methods with square elements (i.e. with df = 0.0) have been presented in [13] and are repeated here for comparison. for coarser meshes in all cases, although optimal convergence is attained at lower refinements as element regularity decreases. In contrast, the new IP methods (Figure 3c) display optimal convergence when ν = 0.49995, irrespective of element regularity, with a small increase in error magnitude as df increases.
The error of the new IP methods is significantly smaller than that of the original IP methods. The SG method with SRI, included here for comparison, also shows optimal convergence, with slightly better accuracy than the IP methods. Finally, the uniform optimal convergence of the new IP methods with respect to λ, independent of df, is illustrated in Figure 3d (shown for the NIPG method). The performance of the other IP methods is the same.

Post-processed stress
For the SG method with Q 1 elements (Figure 5a), for ν = 0.3 there is first-order convergence of the postprocessed stresses, with a slight increase of error magnitude as df increases, while for ν = 0.49995 the convergence rate is poor. For ν = 0.49995, the original IP methods display poor stress convergence rates for coarse meshes; these rates improve with refinement, more quickly for higher df in a given variant ( Figure 5b). In contrast, as for the displacement approximations, the new IP methods produce first-order convergence results uniformly with respect to df for ν = 0.49995 (Figure 5c), with significantly lower error magnitudes that those of the original IP method, and uniformly with respect to the compressibility parameter (NIPG results shown in Figure 5d). The errors of the SG method with SRI are shown in Figure 5c for comparison, and while the convergence rate is first-order, the errors are greater than those of the IP methods, unlike for the displacement approximations.
The contour plots in Figure 6 show, with the exact solution for comparison, the slight deterioration in accuracy of the NIPG post-processed stress field for the component σ xx , for the near-incompressible case, as element regularity decreases. This is vastly improved by mesh refinement. Figures 6g and 6h, when compared to Figures 6e and 6f, illustrate how the contour lines highlight the discrepancies, while the overall stress field is fairly smooth, even for df = 0.3 at mesh 5.  the values of the σ xx field, but appears as 0 when the same scale is used, and at the original scale is improved by refinement (Figure 7d).
The results of the other two IP methods are again similar to those of NIPG.

Displacement approximation
The SG method with Q 1 elements (Figure 8a) shows the same behaviour as for the beam example: optimal convergence when ν = 0.3, with increasing error magnitude for increasing df, and poor convergence, indicating locking, when ν = 0.49995, irrespective of element shape regularity, until high refinement levels. In Figure 8b, optimal convergence is seen for the same elements with SRI applied, for both values of ν and all df, with error magnitudes increasing slightly as df increases. The SG method with Q 2 elements (Figure 8c) shows optimal convergence (second-order, in this case) when ν = 0.3, irrespective of element regularity, though with decreasing accuracy as regularity decreases. With ν = 0.49995, convergence is slower for the coarsest meshes but reaches (or exceeds) optimal convergence quickly with refinement. Results from the original and new IP methods are shown in Figure 9 with performance similar to that in the beam example. In the near-incompressible case, the original IP methods converge poorly for coarse or medium-refinement meshes, depending on the variant, reaching optimal convergence earliest for the least regular meshes (Figure 9a), while the new methods (Figure 9b) converge optimally even at low refinement levels. The error magnitudes are, overall, lower than those of the original methods. SG with SRI is again slightly more accurate in this example. The convergence of the new IP methods with respect to the compressibility parameter is uniform, as shown for NIPG in Figure 9c. The quality of the x-displacement approximation when ν = 0.49995 is shown for SIPG in Figure 10, with a comparison to the exact solution, where it is evident that the accuracy of the method is not lost when general quadrilateral rather than rectangular elements are used. Performance for the component u y is the same, and these results extend to the other two IP methods.

Post-processed stress
The stress L 2 errors for the SG method with Q 1 elements (Figure 11a) show similar convergence behaviour to those in the beam example: with ν = 0.3, convergence is first-order, and for nearincompressibility convergence is poor until high refinements. With SRI (Figure 11b), convergence is firstorder irrespective of the compressibility level, but with ν = 0.49995 the error magnitudes are significantly greater than with ν = 0.3. SG with Q 2 elements (Figure 11c) gives second-order convergence for ν = 0.3; for ν = 0.49995 convergence is also close to second-order except on the coarsest meshes. In most cases, the SG methods show slightly lower accuracy as df increases.
For near-incompressibility, the original IP methods (Figure 12a) give poor convergence for coarse meshes, and the level of refinement at which first-order convergence is attained varies, as in the beam example, with method and element regularity. The new SIPG and IIPG methods (Figure 12d) reach first-order convergence at low refinement levels when ν = 0.49995, while NIPG has rates that vary with element regularity, but with errors smaller than the other two methods because of faster convergence for the coarser meshes. This figure also shows the comparative magnitude for df = 0.3 of the errors of SG Q 1 with SRI and of SG Q 2 , the former significantly larger than the errors of the IP methods and the latter larger until high refinement. Figures 12b and 12c show enlarged plots of the NIPG and IIPG results for in contrast, shows very little deterioration in accuracy, as df increases, the smoothness again increasing with refinement. These different behaviours stem from a lack of smoothness in the field of tr ε for this boundary value problem: although the order of magnitude of tr ε is low, it is amplified by the factor of λ in calculating the direct stresses (see equation (1)). For near-incompressibility, λ is large (O(10 3 )) and thus the contribution of this term (which doesn't appear in the shear component) is significant.   In this problem, with the SIPG method the stabilization parameter values k µ = 10 and k λ = 50 are used. This combination is introduced here for its effect on the stability of the stress results for the new SIPG method.
The initial mesh used is unstructured and graded, and finest around the re-entrant corner (Figure 17a).
While a large portion of the quadrilateral elements in this mesh are close to square, many are of lower regularity and some are very poor in shape (Figure 17b) * . The initial mesh is then refined, and Table 2 details the number of elements and dofs for each method at the various refinement levels. * The measure of deviation in Figure 17b is calculated as 2 divided by the condition number of the weighted Jacobian matrix, by the CUBIT metric "Shape" for quadrilateral elements [10]. Note on domain regularity: The analytical results of [13] for rectangular elements rely on a regularity estimate that assumes a degree of domain regularity not satisfied by the nonconvex polygon of this example (see [5]). (Note that Wihler [23] has developed an extension of this result for general polygonal domains within the framework of weighted Sobolev spaces.) For conforming finite elements, the theoretically predicted rate of convergence for a uniform mesh is 2 3 , based on the interior angle of 3π 2 (see n no. els no. dofs SG Q 1 SG Q 2  IP  1  2303  4796  18802  18424  2  9212  18802  74450  73696  3 36848  74450  296290  294784  4 147392 296290 1182146 1179136  5 589568 1182146 4722562 4716544   Table 2: Mesh details for meshes 1 to 5, for standard Galerkin (SG) and interior penalty (IP) methods for the L-shaped domain example: number of elements and number of degrees of freedom are shown.
[8] and references therein). We, however, use a graded mesh and do not have a theoretical prediction for the convergence rate.
At the origin, there is a stress singularity in the exact solution, which makes this a challenging benchmark boundary value problem for testing computational approximation methods. Figure 18a shows the results of various cases of the SG method. With Q 1 elements, there is a convergence rate of 0.54 for the compressible case. The results are indistinguishable for Q 1 with SRI, and with Q 2 elements the errors are smaller but the convergence rate is the same. For near-incompressibility, the rates decrease significantly for Q 1 elements, and very slightly for SRI. With Q 2 elements, the rate is essentially unaffected but the error magnitudes are higher.

Displacement approximation
The original IP methods (Figure 18b) attain a convergence rate of 0.54, when ν = 0.3, and with ν = 0.49995 the rate is slightly higher although the errors are greater. The new IP methods ( Figure   18c) have a consistent rate of 0.54, for both values of ν, with comparatively little increase in error for the near-incompressible case. Figure 18d shows the higher accuracy of the new methods than of the original for ν = 0.49995, as well as that of the SG method with and without SRI.
The poor displacement approximation of the SG method (x-displacement shown in Figure 19b) when ν = 0.49995 is not an example of locking but of the inaccuracy of the deformation over the domain, and particularly in the region of the origin, as can be seen by comparison to the exact solution in Figure   19a. The three IP approximations in this case are excellent (SIPG result shown in Figure 19c). The corresponding y-displacements display similar accuracy in each case.

Postprocessed stress
With the SG method, convergence for ν = 0.3 is at a rate of about 0.54 for all 3 variants (Figure 20a).
For ν = 0.49995, there is no noticeable decrease in convergence rate (with Q 1 elements there is a slight increase) but the error is several orders of magnitude larger. With the original IP methods (Figure 20b) the convergence rate is likewise 0.54 for ν = 0.3, and there is a decrease in rate for NIPG though not for SIPG and IIPG, but in all three methods there is an increase in error magnitude for ν = 0.49995. Because there is a stress singularity at the origin, an exact stress solution evaluated at nodal values cannot be plotted. Instead, we calculate the exact stress solution at element quadrature points, and find a smooth stress field by projection to the nodal points. This results in the projected exact solution increasing in accuracy with each refinement, most notably with the value at the origin becoming greater with each refinement. Mesh 5 (see Table 2) is used for the contour plots.
The projected exact solution for the σ xx component is shown in Figure 21, with a high stress value appearing at the origin. Figure 22 shows the corresponding SIPG results for various refinement levels. On mesh 1 the results are inaccurate, whether the full domain or the region of the singularity is considered.
Mesh 2 produces an improved approximation across the domain but with lower stresses than the exact solution, and captures the broad features of the exact stress field in the region of the singularity. With the highly refined mesh 5, the full domain is approximated visually accurately, and in the region of the singularity the contour features of the projected exact solution are closely approximated, though small regions of lower stress appear. There are two very small patches of negative stress, above and below the high-stress patch, which do not match the true, all-positive stress field -these patches decrease in area but increase in magnitude with refinement.
The projected exact solution of the σ xy component is shown in Figure 23. The results of the SIPG method are shown in Figure 24, where the behaviour across the full domain is good with both meshes The performance of the SIPG method for the σ yy component (Figure 26) is similar to that for σ xx .
Considering the full domain, by comparison to the projected exact solution in Figure 25a, mesh 1 produces an inaccurate field, mesh 2 gives a significantly improved field and mesh 5 a visually nearlyaccurate one. In the region of the stress singularity, comparison to the projected exact solution in Figure   25b shows that the approximation using mesh 1 is inaccurate, that using mesh 2 gives an improvement, and that using mesh 5 reflects the correct features with high accuracy.
The NIPG and IIPG methods display qualitative behaviour of comparable visual accuracy to the SIPG method in all three components.
The effect of mesh refinement on the postprocessed stress field is marked in this example. However, across the domain broadly, the change in values is small, and the inaccuracies have been highlighted in presentation of the results by the choice of contour values to enhance understanding of the behaviour of the methods. In the immediate region of the origin, the maximum and minimum values of stress are produced, and the highly magnified images show the mediocre performance of the IP methods in attaining these extremes using mesh 1, and the vastly improved performance using more refined meshes. The exact solution is given in [13].
In this problem, with the SIPG method the stabilization parameter values k µ = 50 and k λ = 250 are used. These values are introduced here for their effect on the stability of the displacement results for the original and new SIPG methods.    Stabilization parameters higher than the usual choices were required for stability of the SIPG method with ν = 0.49995, in this example. While the question of how to choose stabilization parameters is one that is relevant to IP methods broadly and is an area of study in its own right, we note in the context of this work that increasingly high values of k λ are required for stability here as element shape regularity decreases, i.e. as df increases, as demonstrated in Figure 28.
A contour plot of an oblique cross-section of the cube shows by comparison with the exact solution that the accuracy of the x-displacement approximation of the SIPG method with cubic elements (df = 0.0) is maintained for a mesh of distorted elements (with df = 0.3) when ν = 0.49995 (Figure 29). The accuracy of the y-displacement approximation is similar, and these results extend to the other two IP methods.

Postprocessed stress
As in the first two examples, the L 2 error of the postprocessed stress of the SG method with Q 1 elements shows first-order convergence for all df when ν = 0.3, with a slight decline for df = 0.3, and diverges when ν = 0.49995 (Figure 30a). When SRI is applied, first-order convergence is displayed for all df for both values of ν, with larger errors for the nearly incompressible case (Figure 30b). The original IP methods show diverging stress L 2 errors for the nearly incompressible case (Figure 30c). All three new IP methods tend towards first-order convergence for all df when ν = 0.3; when ν = 0.49995 the rates are lower than first-order, declining with increasing df and poorest for SIPG, but improve with refinement (Figures 30d, 30e and 30f). The error magnitudes are larger for the nearly incompressible case than for the compressible case, with decreasing element regularity also resulting in greater errors.
The IP errors are, however, smaller than those of the SG with SRI.
The qualitative behaviour of the SIPG method on mesh 5 for the individual components σ xx and σ xz , for ν = 0.49995, is shown in Figures 32 and 33, with the exact solutions shown in Figure 31 for comparison.
There is a clear deterioration in accuracy in the approximate σ xx field as df increases: with df = 0.0, the result is visually accurate, with df = 0.1 there is a loss of smoothness of the field, and with df = 0.3 the field is significantly more patchy, although the general areas of high and low stresses are still reflected. The coarser the meshes, the less clearly these zones are displayed (images not shown here), indicating that refinement improves the quality of the postprocessed field in this component. The results for the direct stresses σ yy and σ zz are similar. The approximate shear stress component σ xz (Figure 33) shows comparatively little deterioration in accuracy, though some loss of smoothness, as df increases, as do the shear components σ xy and σ yz .
The other two IP methods similarly display increasing "patchiness" in the direct stresses as df increases, and similarly maintain good overall approximations, though with loss of smoothness as df increases, in the shear components.

Conclusion
The first and primary aim of this paper was to establish that the uniform convergence, with respect to the compressibility parameter, of the new IP methods of [13] extends to the use of general quadrilateral/hexahedral elements. In all four model problems presented here, this uniform convergence is indeed seen. While the IP approximation errors are slightly larger in many cases for greater degrees of distortion, this is not a feature of near-incompressibility only: this behaviour is seen in the case ν = 0.3 as well, and is also displayed by the SG method for ν = 0.3. Contour plots of the displacement approxi-mations indicate that the use of general quadrilateral/hexahedral, instead of rectangular, elements does not visibly diminish the quality of the displacement results.
The second aim was to ascertain the extent to which a stress field postprocessed from the IP displacement approximation is accurate, particularly in the near-incompressible case and for non-rectangular elements.
We see in the examples of the cantilever beam and the square plate that, for the stress approximation, the first-order convergence rate of the L 2 error achieved for the case of ν = 0.3 using all three IP methods is also achieved for the case of ν = 0.49995, irrespective of element shape, except in the case of the NIPG method for the square plate, where the error is nevertheless smaller than for the other IP methods. In the example of the L-shaped domain, the convergence rate is lower than first-order for ν = 0.3 and this rate is exceeded for ν = 0.49995. For the example of the cube, convergence is first-order when ν = 0.3, and for ν = 0.49995, convergence tends towards or exceeds this rate.
These examples strongly suggest that the stress error of the new IP methods converges uniformly in the L 2 -norm with respect to the compressibility parameter.
To assess the quality of the stress approximation on a given mesh (i.e. at a given refinement level), a continuous field was obtained for visualisation by a projection of the stress values calculated at quadrature points onto the mesh nodes, and compared component-wise to the analytical stress expression through contour plots. In general, with ν = 0.49995 and general quadrilateral/hexahedral elements, IP results in the direct stress components were of lower quality than with rectangular elements, and more evidently so the lower the element shape regularity, while in the shear stress components the approximate fields lost only smoothness. With sufficient refinement, the quality was recovered, but the refinement levels necessary for good-quality stress approximations corresponded in some cases to very large systems of equations. The displacement approximations typically give visually accurate contour plots at several refinement levels lower than the postprocessed stresses do, for a given problem.
This phenomenon is not, however, unique to the IP methods. For example, the SG method with Q 2 elements likewise produces low-accuracy stress fields in the near-incompressible case at a refinement level at which the displacements are extremely good, as the error plots show. Poor performance in the stress in the near-incompressible case is, moreover, not limited to the case of non-rectangular elements.
Obtaining accurate results from postprocessing would thus seem to require alternative strategies, even where error convergence rates are very good.
A significant addition and useful complement to the computational investigation presented here would be a theoretical displacement error analysis of the new IP formulation for the case of general quadrilateral elements. More generally, there is scope for an extension of this formulation and the corresponding analyses to nonlinear problems. Regarding a numerical, implementational aspect of this work, direct solvers were used in the examples in this study, but are computationally inefficient for large systems.
An investigation into the design of an appropriate preconditioner to be used with an iterative solver in the near-incompressible regime, with the aim of decreasing computing time and memory usage, would therefore be of value.