1 Introduction

Super-convergence computation in the finite element method (FEM) has been a research hotspot for decades. Super-convergent solutions usually gain higher accuracy compared with ordinary finite element (FE) solutions on the same mesh, and can be used to construct a posteriori error estimator in adaptive FE analyses. One type of a posteriori error estimate is residual-based^[1], which estimates errors on each element using residuals of the differential equations. Another type is recovery-based, which smooths the gradients of FE solutions across element boundaries, including the well-known super-convergent patch recovery (SPR) method^[2-3] and some other related methods^[4-6]. Most of the above methods estimate errors in L₂ norm rather than the more intuitive and stronger maximum norm, which requires point-wise super-convergent solutions.

The element energy projection (EEP)^[7] is a technique for point-wise super-convergence computation in the post-processing stage of FEM, which assumes the well-known energy projection theorem^[8] to almost hold true on a single element of degree m (m > 1). This assumption proves to be true for the one-dimensional (1D) FEM, and the consequent EEP technique has been applied to various 1D problems^[9-10]. Moreover, it has been proved mathematically for the 1D FEM that the EEP solutions of simplified form achieve convergence at least one order higher than their FE solutions. However, in the case of two-dimensional (2D) or three-dimensional (3D) problems, the assumption behind the EEP technique is no longer valid, and hence the EEP formulae cannot be derived directly. Thus, whether the 1D-based EEP technique can be somehow extended and applied to multi-dimensional problems remains an open question and a severe challenge. The present paper attempts to give a positive answer to the open question in an abstract and conceptual way.

An intuitive idea to conquer this challenge is to treat multi-dimensional problems as generalized 1D problems, so that the combination of 1D FE discretization with the 1D-based EEP technique can be applied. This idea is analogous to the calculation of a triple integral,

(1)

where the volume integral is recursively calculated via three 1D integrals, i.e., integrating in a dimension-by-dimension (D-by-D) manner. Similarly, the key to achieve super-convergence computation of 3D problems is to implement EEP formulae in a recursive D-by-D manner, which calls for a unified EEP formulation for each generalized 1D problem.

This paper starts by introducing the so-called generalized 1D problem, generated by discretizing only in one dimension of a multi-dimensional problem with the generalized 1D FE technique. This allows a unified 1D-based EEP formulation to be developed, and a unified derivation of the unified EEP formulae is given. In the generalized 1D FEM and EEP formulae, unknowns on generalized nodes are discretized recursively using the generalized 1D FEM, until a fully discretized FE model is reached, which is equivalent to that from a full FE discretization. The recursive D-by-D discretization leads to a reverse D-by-D recovery procedure, which is the major gain of this paper. A critical issue in the recovery procedure is how to extract the generalized nodal solutions from the FE solutions, about which, by intuitive and numerical analyses, a set of workable replacing solutions have been found out and are suggested explicitly, which can render super-convergent solutions of displacements and first derivatives over the 3D domain. Numerical examples are given to verify the feasibility and effectiveness of the strategy proposed in this paper.

2 Generalized 1D problem 2.1 Model problems

Consider the following 1D second-order linear ordinary differential equation (ODE) and the multi-dimensional Poisson's equation as a batch of the model problems:

(2a)

(2b)

(2c)

with appropriate boundary conditions, where p (p > 0) and q are functions of x, and Ω is a 2D or 3D domain. Here and below, and so are other subscripts. L is the associated differential operator. Function u will be called displacement, and f is the load term. The corresponding energy inner products of model problems (2) are

(3a)

(3b)

(3c)

The batch model problems in Eq. (2) are to be discretized into generalized 1D problems by the generalized 1D FEM in the following.

2.2 Generalized 1D FEM

Generalized 1D FEM discretizes multi-dimensional problems in only one dimension by FE approximation and leaves the continuous unknowns of other dimensions defined on generalized nodes of the 1D FE model. Thus, the unknowns for 1D problems are simply nodal values, and those for 2D and 3D problems are extended to be functions on nodal lines and nodal faces respectively. For convenience, these three types of generalized 1D methods are named FEM, FEM of lines (FEMOL)^[11], and FEM of faces (FEMOF), respectively, and their generalized 1D element models in local coordinates are demonstrated in Fig. 1.

For 1D problems, the generalized 1D FEM is simply the conventional 1D FEM as shown in Fig. 1(a), in which the trial function on an element can be written as

(4)

where [N(ξ)] is the element shape function matrix, and {d^N}^e is the element nodal displacement vector. Eventually, a global stiffness equation of the following form can be derived:

(5)

where [K] and {P} are the global stiffness matrix and the load vector, respectively. Equation (5) is a linear algebraic equation and can be solved as in the conventional FEM.

For 2D problems, the generalized 1D FEM has been named FEMOL^[11] as mentioned above. Figure 1(b) shows a typical FEMOL element, and the trial function can be written as

(6)

where {d^L(η)}^e is the element nodal line displacement vector. Using such trial functions and applying the minimal potential energy theorem (or virtual work equation), a system of ordinary differential equations (ODEs) called FEMOL ODEs of the following form can be derived:

(7)

For the associated boundary conditions and the related notations, we can refer to Ref. [11]. It has been mathematically proved^[12] that the FEMOL solution obtained via the solution of Eq. (7) has a convergence order of h^m+1, where h is a measure of the largest size of the FEMOL elements in the discrete direction.

For 3D problems, the generalized 1D FEM is conceptually named FEMOF. Figure 1(c) shows a typical FEMOF element, and the trial function can be written as

(8)

where {d^F}^e is the element nodal face displacement vector. By using the same means as in 2D problems, a system of partial differential equations (PDEs) called FEMOF PDEs of the following expanded form can be derived:

(9)

with associated boundary conditions, where the global nodal face displacement {d^F} is formed by assemblage of {d^F}^e, and [B_i] (i=1, 2, …, 6) and {F^F} are functions of η and ζ, respectively. Note that Eq. (9) is a system of PDEs defined on nodal faces, and with its solutions being somehow obtained, Eq. (8) gives the FEMOF solution for the 3D problem in Eq. (2c).

The features of the three types of generalized 1D FEM are summarized in Table 1.

3 EEP formulation 3.1 Unified EEP formulae of generalized 1D FEM

The generalized 1D FEM for problems of different dimensions allows the 1D based-EEP technique to be adopted to derive a set of unified EEP formulae on a typical element, which can be cast into the following form:

(10a)

(10b)

where (·)^* indicates the EEP super-convergent solution, ξ_a∈[ξ₁, ξ₂]≡[-1, 1], N₁ and N₂ are two linear shape functions associated with the two generalized end nodes, J is the determinant of the Jacobian of element geometric mapping, and other specific notations are described in Table 2. It is seen that both super-convergent solutions are corresponding FE solutions modified by adding two integral recovery terms. The residual term (f-Lu^h) in these integrals can be calculated when the solutions of Eq. (5), Eq. (7), or Eq. (9) are obtainable.

For 1D problem (2a), the EEP formulae (10) for super-convergent displacement u_a^* and first derivative (u_n^*)_a≡u_a'^* at an interior point ξ_a of an element of length h are equivalent to those in Ref. [7], except for different coordinate systems. The residual term (f-Lu^h) involves nodal values {d^N} and can be calculated directly. Also, it has been proved mathematically that for m > 1, the displacement u_a^* and derivative u_a'^* in Eq. (10) have super-convergent orders of at least h^m+2 and h^m+1, respectively, which are both one order higher than their corresponding FE solutions.

For 2D Poisson's equation (2b), the EEP formulae (10) are identical to those in Ref. [13] with the expression for the ξ-normal derivative being

(11)

where J_L, c, and b are geometric mapping terms of the related FEMOL element^[12]. Here, the residual term (f-Lu^h) involves {d^L}, {d^L'}, and {d^L"}, which can be obtained in conventional FEMOL by solving the FEMOL ODEs of Eq. (7) using existing ODE solvers. Numerical examples have shown that for m > 1, and calculated by Eq. (10) have convergence orders of h^m+2[13] and h^m+1, respectively, which are both one order higher than their corresponding FEMOL solutions.

For 3D Poisson's equation (2c), the EEP formula (10a) for displacement is outlined in Ref. [14] and formula (10b) can be calculated similarly, with the expression for ξ-normal derivative being

(12)

where J_F, a₁, a₂, and a₃ are geometry-related quantities like those in Eq. (11). The residual term (f-Lu^h) involves {d^F} and its various derivatives {d_η^F}, {d_ζ^F}, {d_ηη^F}, {d_ηζ^F}, and {d_ζζ^F}. If the true solutions of Eq. (9) could be obtained somehow, Eq. (10) would also be expected to give super-convergent solutions of one order higher than their corresponding FEMOF solutions.

3.2 Unified derivation of EEP formulae

This subsection gives a unified derivation of Eq. (10) based on 1D EEP technique.

The linear (in the ξ-direction) test functions are taken as the following forms:

(13a)

(13b)

(13c)

In the following derivation, the 3D case is taken as the working example with and replaced by u^h and v^h for short. The derivations in 1D and 2D cases can be done in an almost identical way with the following replacements:

(14a)

(14b)

By substituting the trial function u^h and the test function v^h of the FEMOF into the energy inner product in Eq. (3), the EEP becomes

(15)

where e^h=u-u^h, and Ω_e is the domain of a typical FEMOF element with F₁ and F₂ being its two end nodal faces, respectively. Also, divide Ω_e into two subdomains Ω₁ and Ω₂ by the interface F_a defined by ξ =ξ_a.

Firstly, by taking integration by parts in different manners on each of the two subdomains, a^e(e^h, v^h) can be written as

(16)

Since v^h on the four edges (η_j = ±1, ζ_j = ±1) is independent of its values on the interior domain defined by [ξ₁, ξ₂]×(η₁, η₂) × (ζ₁, ζ₂), the end side terms in Eq. (16) can be omitted which lead to the following expanded form:

(17)

Secondly, noticing that a^e(e^h, v^h) is the summation of two terms, a^e(e^h, N₁v₁^h) and a^e(e^h, N₂v₂^h), and taking integration by parts in different manners for each of the two terms, lead to another expanded form,

(18)

where the end side terms have been omitted again. Equating Eq. (17) and Eq. (18) yields the following relationship:

(19)

Some trivial rearrangement and reformation of Eq. (19) give

(20)

with

(21)

where ε is the error term and is conjectured to be negligible.

On one hand, by choosing a particular test function v^h =v₀(η, ζ)·(ξ-ξ_a), which implies , , and ignoring ε term, Eq. (20) becomes

(22)

where dV=Jdξdηdζ, dA=J_Fdηdζ, and Le^h=f-Lu^h. From the arbitrariness of v₀ (η, ζ) on (η₁, η₂) × (ζ₁, ζ₂), Eq. (10a) in the 3D case is arrived at.

On the other hand, by choosing another particular test function v^h =v₀ (η, ζ), which implies v_a^h = v₁^h = v₂^h = v₀ and vn^h|_a = ((a₂v_{0, η} + a₃v_{0, ζ})/J_F)_a, Eq. (20) becomes

(23)

where the end side terms have been omitted again. Ignoring ε term and from the arbitrariness of v₀(η, ζ) on (η₁, η₂) × (ζ₁, ζ₂), the following expression is derived:

(24)

The last four terms (e_η^ha₂ + e^ha_{2, η}+e_ζ^ha₃+e^ha_{3, ζ})_Fa in Eq. (24) are conjectured to be of higher convergence orders and are hence neglected, which leads to Eq. (10b) in the 3D case.

Note that the EEP formulae of different dimensions are not limited to a single equation like Eq. (2). In fact, EEP formulae for a system of ODEs have also been developed with satisfactory results^[17], and those for a system of 2D PDEs can also be derived in a similar procedure.

4 Super-convergence of multi-dimensional FEM 4.1 Equivalent discretization and recovery

The generalized 1D FEM and the unified EEP formulae for 2D or 3D problems leave a system of FEMOL ODEs or FEMOF PDEs to be solved on generalized nodes. In this paper, those ODEs or PDEs will be further discretized using the generalized 1D FEM again and again recursively, until ultimately a fully discretized standard FE model is achieved. Along with each 1D discretization, the corresponding EEP formulae, derived in an analogous way with that in Subsection 3.2, can be applied. Note that this D-by-D discretization is conceptual and is not practically implemented since the final FE model can be directly achieved by full FE discretization. However, as a reverse procedure of the D-by-D discretization, a D-by-D recovery strategy for recursive super-convergence computation is readily obtainable and developed.

Taking 3D problem as an example, the upper part of Fig. 2 shows the equivalent discretization of 3D FEM in a D-by-D manner. The meshes conceptually generated by three successive 1D discretization steps can be equivalently achieved by a direct full discretization, with the only requirement being that the 3D domain is partitioned by transverse faces in three directions as shown in the rightmost of Fig. 2. The 3D FE meshes thus generated are called quasi-FEMOF meshes. Thus, the trial function on a 3D element is naturally defined as

(25)

which is identical to that derived by three successive steps of 1D discretization. Similarly, the equivalent discretization also leads to identical global stiffness equations.

4.2 Extraction of generalized nodal solutions

The lower part of Fig. 2 sketches the flowchart of the D-by-D recovery procedure, which is basically recursive application of the unified EEP formulae. The key issue in its implementation is how to deal with the generalized nodal solutions, including displacement and various derivatives. The true solutions on the generalized nodes are assumed to be available in those EEP formulae, meaning that FEMOL ODEs or FEMOF PDEs have to be solved exactly, which is beyond the intention of this paper. With FE solutions at hand and various EEP formulae available to help extract super-convergent solutions at generalized nodes, a natural and logical idea is to use certain kind of combination of FE and EEP solutions to replace the true solutions of systems of ODEs or PDEs. With intuitive analyses and numerical experiments, this paper has found out, among others, an appropriate set of replacing solutions used as the generalized nodal solutions, which can ideally render overall super-convergence of one order higher than FE solutions.

Table 3 presents the proposed set of replacing EEP solutions, and Table 4 and Table 5 give further explanations for 2D and 3D problems, respectively.

For 2D problems, the above EEP-based super-convergent computation has been around for some time and has been successfully applied to self-adaptive analysis of 2D FEM^[18]. But for a long time, the success for 2D problems has not been able to be extended to 3D problems, especially on irregular meshes, since the extraction of nodal face solutions involves far more terms and choices. The choice shown in Table 5 has been found to work well for 3D problems and makes the 1D-based EEP technique applicable in a unified fashion from 1D to 3D problems throughout. Further for 3D problems, apart from the super-convergent displacement and normal derivatives in Eq. (10), the other two super-convergent first derivatives can also be obtained by interpolating along the ξ-direction to {d_ζ^F}^* and {d_η^F}^* on nodal faces. Using these terms, super-convergent ξ derivatives can be derived from Eq. (12). Such implementation is applicable for not only EEP formulae of 3D Poisson's equation but also for those of other 3D linear problems, and has the potential to work for EEP formulae of even higher dimensional problems.

5 Numerical examples

In this section, 3D problems including Poisson's equation and elasticity problem are both taken as numerical examples to demonstrate the feasibility and effectiveness of the proposed strategy. Examples of 1D and 2D problems will not be given because their super-convergent solutions are the prerequisites for the computation of 3D problems.

The main purpose to conduct numerical tests is to check errors in various solutions based on FE results on quasi-FEMOF meshes. In doing so, the maximum norm is adopted to calculate the errors. In general, these errors can be cast into the form of Ch^ρ, where h is the largest element size, and C is a constant independent of h. Given the errors on a series of uniformly refined meshes (i.e., h is halved for all elements in each refinement), the convergence order ρ can be approximately estimated by with r being the ratio of the two errors before and after refinement. In addition, let ρ denote the asymptotically converged value of ρ and N_e denote the number of elements used.

5.1 3D Poisson's equation with homogeneous boundary conditions

Consider the following 3D Poisson's equation:

(26)

which is defined on an irregular hexahedron domain Ω as shown in Fig. 3. The load term f is determined such that the true solution u is

(27)

Elements of degree m(m=2, 3, 4) are used in computation, and linear elements have been excluded because no super-convergence orders can be gained. The following errors have been calculated:

(28)

where u^h and u_ξ^h, u_η^h, u_ζ^h denote 3D FE solutions of displacement and its derivatives, and u^* and u_ξ^*, u_η^*, u_ζ^* denote the related EEP solutions.

By computing on a series of uniformly refined meshes, like the one shown in Fig. 4, the convergence orders for displacements and derivatives in three directions are given in Table 6 and Table 7.

5.2 3D elasticity problem with inhomogeneous boundary conditions

In a typical 3D elasticity problem, the basic unknown is the displacement vector {u} = {u v w}^T, and the displacement-based equation has the form

(29)

where [S] is a 6 × 3 linear differential operator matrix, [D] is a 6 × 6 stiffness matrix, and {f} is the corresponding body force vector.

Consider the 3D elasticity problem with the equation in Eq. (29) defined on domain Ω, which is the part of a unit cubic above the plane x + y = 10z, as shown in Fig. 5.

The load vector {f} is determined such that the true solution is

(30)

The boundary conditions are set to be

(31)

where {q}≡[n][D][S]{u} is the surface traction of an arbitrary plane with {n} being its 3 × 6 out normal direction matrix, and {u} is defined in Eq. (30).

Elements of degree m(m= 2, 3) are used in computation, and the following errors of displacements and ξ-derivatives have been calculated:

(32)

and the errors e_η^h, e_η^*, e_ζ^h, and e_ζ^* are calculated in a similar manner.

By computing on a series of uniformly refined meshes, the convergence orders for displacements and derivatives in three directions are given in Table 8 and Table 9.

It is noted that when m=2 the three derivatives calculated with EEP technique do not seem to be super-convergent, although their errors are siginificantly smaller than those of FE results. One of the discernable reasons is the difference between the true surface traction {q} on the inhomogeneous boundary and its interpolation over discrete FE nodes. Despite the difference in the vicinity of the boundary (i.e., the so-called boundary layers), the derivatives calculated with EEP technique in the majority of the domain are super-convergent. Let Ω^* be a subdomain of Ω excluding 1/8 volume of Ω near the boundary x + y = 10z. By re-calculating errors e_ξ^*, e_η^*, and e_ζ^* according to Eq. (32) with Ω substituted by Ω^*, convincing results of super-convergent derivatives are achieved as shown in Table 10.

From the above numerical examples and all the examples not given here, the convergence orders of displacement and derivatives can be summarized as shown in Table 11. In summary, the proposed EEP technique for 3D FEM gains super-convergence at least one order higher than the ordinary FE solutions.

6 Concluding remarks

In this paper, multi-dimensional problems in FEM are viewed as a result of recursive generalized 1D FE discretization, so that the 1D-based EEP recovery technique can be effectively used in a recursive manner. With proper extraction of the generalized nodal solutions required in EEP formulae, point-wise super-convergent solutions of displacements and three first derivatives are obtained for 3D FEM on quasi-FEMOF meshes. Numerical examples have verified the feasibility and effectiveness of the proposed strategy. Recent studies have also shown that the super-convergent solutions are well applicable in adaptive analysis of 3D FEM using maximum norm, and the latest progress will be presented in other publications.