Non-uniform order mixed FEM approximation - Semantic Scholar
Journal of Computational Physics 231 (2012) 436–453
Contents lists available at SciVerse ScienceDirect
Journal of Computational Physics journal homepage: www.elsevier.com/locate/jcp
Non-uniform order mixed FEM approximation: Implementation, post-processing, computable error bound and adaptivity Mark Ainsworth ⇑,1, Xinhui Ma 2 Mathematics Department, Strathclyde University, 26 Richmond Street, Glasgow G1 1XH, Scotland, United Kingdom
a r t i c l e
i n f o
Article history: Received 29 March 2011 Received in revised form 31 August 2011 Accepted 8 September 2011 Available online 21 September 2011 Keywords: Posteriori error estimation Mixed finite element method Computable error bounds Electromagnetics Flow in porous media Post-processing
a b s t r a c t The present work provides a straightforward and focused set of tools and corresponding theoretical support for the implementation of an adaptive high order finite element code with guaranteed error control for the approximation of elliptic problems in mixed form. The work contains: details of the discretisation using non-uniform order mixed finite elements of arbitrarily high order; a new local post-processing scheme for the primary variable; the use of the post-processing scheme in the derivation of new, fully computable bounds for the error in the flux variable; and, an hp-adaptive refinement strategy based on the a posteriori error estimator. Numerical examples are presented illustrating the results obtained when the procedure is applied to a challenging problem involving a ten-pole electric motor with singularities arising from both geometric features and discontinuities in material properties. The procedure is shown to be capable of producing high accuracy numerical approximations with relatively modest numbers of unknowns. Ó 2011 Elsevier Inc. All rights reserved.
1. Introduction Adaptive high order finite element schemes incorporating both local mesh refinement and variable order approximation have proved to be one of the most effective methods for the numerical treatment of problems governed by elliptic partial differential equations [1–3], and are typically referred to as the hp-version of the finite element method. The solutions of such problems typically exhibit local singularities in the neighbourhood of re-entrant corners, points at which there is a change in type of boundary condition, and at interfaces and cross-points where the material properties undergo a jump change in value. The theory and implementation of hp-FEM is reasonably well-understood and widely practised in the case of H1-conforming finite element approximation. However, many practitioners prefer to use a mixed finite element scheme [4] for the approximation of wide classes of problems arising in fluid flow, electromagnetics and elastostatics where one is often more concerned with the flux or stress than in the potential or displacement. Although hp-FEM has been developed and analysed for such problems and shown to produce exponential rates of convergence for problems with the same kinds of singularities described above [5–7], the practical uptake of such methods is comparatively low in comparison with the interest seen for H1-conforming approximation. In search of reasons for the relatively poor uptake, one may look to the fact that the underlying function spaces are less straightforward to discretise in terms of: the construction of appropriate families of hierarchical basis functions; the assembly of the stiffness matrices
⇑ Corresponding author. 1 2
E-mail address: [email protected] (M. Ainsworth). The work of MA was supported in part by EPSRC under grant ‘‘Numerical Algorithms and Intelligent Software for the Evolving HPC Platform’’. The work of XM was supported by the EPSRC grant (EP/E040993/1) ‘‘Adaptive Numerical Methods for Optoelectronic Devices’’.
0021-9991/$ - see front matter Ó 2011 Elsevier Inc. All rights reserved. doi:10.1016/j.jcp.2011.09.011
M. Ainsworth, X. Ma / Journal of Computational Physics 231 (2012) 436–453
437
and load vectors in the case of non-uniform local orders of approximation; and, the implementation of a fully reliable a posteriori error control and adaptive strategy is more complicated. The present work aims to provide the practitioner with a straightforward and focused set of tools and corresponding theoretical support for handling each of these issues. We begin with formulating a problem in mixed form and describing the details of the discretisation of the problem using non-uniform order mixed finite elements of arbitrarily high order. Full details of the necessary procedure are supplied along with representative examples. We then turn to the issue of post-processing of non-uniform order mixed finite element approximations. Post-processing schemes for low order mixed finite element schemes of uniform order have been investigated by a number of authors [8,9]. Here, we generalise these schemes to the case of high order non-uniform order approximation. The schemes produce a smooth (H1-conforming) approximation to the primal unknown through the successive solution of an elementwise local Neumann problem, followed by a local inter-element smoothing and finally, the solution of an elementwise local Dirichlet problem. The local nature of the post-processing means that the cost is negligible in comparison with that entailed in the computation of the mixed finite element approximation and, in addition, is readily implemented. This post-processing procedure is then exploited in the construction of an a posteriori estimator for the error in the flux variable. Although a number of a posteriori error estimators have been proposed for mixed finite element approximation [10–15] none of these estimators produce computable guaranteed bounds for the error in the case of non-uniform order mixed finite element approximation. We extend the analysis presented in [16] for the case of lowest order mixed finite element approximation to the case of non-uniform order approximation and derive a fully computable bound for the error in the flux variable. The cost of computing the error bound is negligible once the post-processed solution is in hand. Finally, we describe an hp-adaptive refinement strategy based on the a posteriori error estimator. The present paper contains all of the details and theory needed to implement a non-uniform order adaptive mixed finite element scheme for the solution of problems in mixed form along with guaranteed error control. Numerical examples are presented illustrating the results obtained when the procedure is applied to the solution of a simple L-shaped model problem, and for a challenging problem involving a ten-pole electric motor with a number of singularities arising from both geometric features and discontinuities in material properties. The procedure is shown to be capable of producing high accuracy numerical approximations with relatively modest numbers of unknowns. 2. Model problem Consider the following model problem:
r A grad u ¼ 0; div r þ f ¼ 0 u ¼ uD nr¼g
in X ð1Þ
on CD on CN ;
where X is a bounded Lipschitz domain in the plane with outer unit normal n on the polygonal boundary C = CD [ CN split into a relatively open Neumann boundary CN and a closed Dirichlet boundary CD. The data f 2 L2(X), g 2 L2(CN) and uD 2 CðCD Þ are given. We assume uD is piecewise polynomial for simplicity. The space L2(X) consists of functions which are square integrable, while the space H1(X) in addition requires all components of the gradients are square integrable. The matrix-valued function A 2 L1 ðX; R22 Þ is assumed to be symmetric positive definite and, for ease of exposition, is assumed piecewise constant on sub-domains of X. The mixed variational formulation of the problem can be written as: find (u, r) 2 L2(X) H(div; X) such that n r = g on CN and
ðA1 r; sÞ þ ðu; div sÞ ¼
R CD
uD n s ds;
ðdiv r; v Þ þ ðf ; v Þ ¼ 0; R where ðv ; sÞ 2 L2 ðXÞ HD ðdiv; XÞ; ðf ; v Þ ¼ X f v . The space H(div; X) denotes the space of vector fields {s 2 L2(X)2 : div s 2 L2(X)}, which we equip with the weighted inner product defined by
ðr; sÞHðdiv;XÞ ¼ ðr; sÞA1 þ ðdiv r; div sÞ; where ðr; sÞA1 ¼ ðA1 r; sÞ, with the associated norms denoted by k kH(div;X) and k kA1 . When x 2 X we shall use the notation (, )x to denote the integral inner product over x and omit the subscript in the case where x is the physical domain X. The subspace HD(div; X) is defined to be {s 2 H(div; X) : n s = 0 on CN}. For discretization, we consider a family of partitions P of the domain X into the union of non-overlapping, shape-regular triangular elements such that the nonempty intersection of a distinct pair of elements is a single common node or single common edge, and such that CD \ CN is a subset of the nodes in P. Each element is assumed to be shape regular, but the mesh maybe highly refined locally. Each element K is assigned a local polynomial degree pK. The variational problem can be discretized by introducing finite dimensional subspaces M P L2 ðXÞ and Q P Hðdiv; XÞ relative to the partitions P as follows:
438
M. Ainsworth, X. Ma / Journal of Computational Physics 231 (2012) 436–453
MP ¼
v 2 L2 ðXÞ : v jK 2 Pp
K1
8K 2 P ;
and
Q P ¼ s 2 Hðdiv; XÞ : sjK 2 PpK PpK 8K 2 P ; where Pk denotes polynomials of total degree at most k. The space Q P corresponds to the Brezzi–Douglas–Marini finite element [4] in the case pK > 0, whilst the lowest order Raviart–Thomas finite element is used in the case pK = 0. The pair of spaces MP and Q P provides a stable approximation of the mixed problem (1). FE We seek a pair uFE 2 M P Q P satisfying n rFE P ; rP P ¼ PN g on CN
R FE A1 rFE P ; s þ uP ; divs ¼ CD uD n s ds; div rFE P ; v þ ðf ; v Þ ¼ 0 for all ðv ; sÞ 2 M P Q P \ HD ðdiv; XÞ. Here, PNg is the piecewise orthogonal projection onto the space PpK ðcÞ for each edge c CN \ @K and all K 2 P. We introduce hierarchic bases M P ¼ spanfwm : m 2 Mg L2 ðXÞ and Q P \ HD ðdiv; XÞ ¼ spanfun : n 2 Ng HD ðdiv; XÞ, which will be defined in next section. The corresponding approximations to r and u are given respectively by
rFE P ¼ uFE P ¼
P n2N
rn un þ rNP ;
P
m2M
um wm :
where rNP 2 Q P is any function satisfying n rNP ¼ PN g on CN . The coefficients are determined by the solutions of the linear system
P n2N
P
n2N
rn A1 un ; ui þ
P
um ðwm ; divui Þ ¼
m2M
rn ðdivun ; wj Þ ¼
R
X
f wj
R CD
uD n ui A1 rNP ; s 8i 2 N;
8j 2 M;
or in matrix form
Br þ Cu ¼ b >
C r¼f
;
ð2Þ
where, for i, n 2 N and j, m 2 M
Bin ¼ ðA1 un ; ui Þ; C im ¼ ðwm ; divui Þ; R bi ¼ CD uD n ui A1 rNP ; ui ; R fj ¼ X f wj : 3. Basis functions and system assembly in the case of locally non-uniform order approximation In this section, we discuss the details of the calculation needed to assemble the matrices appearing in (2) 3.1. Hierarchic basis The reference domain is taken to be the right triangle
b ¼ fð^x; y ^Þ : 1 < y ^ < 1; 1 < ^x < y ^g K with vertices and edges labeled as shown in Fig. 1. Let ^ k0 ; ^ k1 and ^ k2 denote the usual barycentric, or area coordinates, given by
^k0 ¼ 1 ð^x þ y ^Þ; 2
^k1 ¼ 1 ð^x þ 1Þ; 2
^k2 ¼ 1 ðy ^ þ 1Þ: 2
Each edge is described by a pair of vertex numbers ½ i may be assigned a unique parametrization given by
^nij ¼ ^kj ^ki 2 ½1; 1 ij 2 f01; 12; 20g:
j and oriented from vertex i to j in anti-clockwise sense. An edge
M. Ainsworth, X. Ma / Journal of Computational Physics 231 (2012) 436–453
439
Fig. 1. Notation for reference triangle element and orientation of edges.
For later reference, we introduce bubble functions on each edge and triangle by
bij ¼ ^ki ^kj
on edge ij 2 f01; 12; 20g
and
b012 ¼ ^k0 ^k1 ^k2 : Note that these functions vanish on the boundary of the entity with which they are associated. A hierarchic basis on a reference interval I = [1, 1] can be constructed by a set whose lth member function is a polynomial of degree l = 0, 1, . . .. This set can be used to construct general bases on a triangle and is referred to as the primary basis. One obvious choice for the functions would be the monomials of degree l. However, in practice, the function is taken to be the Legendre polynomial Ll of degree l, defined inductively by the relations
L0 ðsÞ ¼ 1;
L1 ðsÞ ¼ s
and
Llþ1 ðsÞ ¼
2l þ 1 l sLl ðsÞ Ll1 ðsÞ; lþ1 lþ1
l ¼ 1; 2; . . .
The primary basis for an interval I may be used to construct polynomial bases on edges by taking the arguments of the primary basis functions to be edge parametrization nij 2 [1, 1]. The functions defined above are used to build discrete subspaces of L2, H1 and H(div) as follows. 3.1.1. L2 space The space L2 imposes no continuity across element interface. This dramatically simplifies the construction of finite eleb is defined by ^ I of order p on the reference triangle K ment space, and a hierarchic basis functions w l;m
^ I ¼ Ll ð^n20 ÞLm ð^n01 Þ 0 6 l; m; l þ m 6 p: w l;m
ð3Þ
The set of hierarchic basis functions consists entirely of interior functions. The term interior reflects the nature of the degrees of freedom associated with the functions. The absence of inter element continuity requirements for the L2 space means the coefficients of the functions can be chosen independently on all elements, i.e., all degrees of freedom are internal to the element. To transform the local basis on the reference element to the global basis on a physical element, let the mapping from the b to the physical element K be denoted by FK : K b ! K. We assume that FK is a differentiable bijection, and reference element K let JK = DFK denote the Jacobian matrix of the transformation. b Þ defined on ^ 2 L2 ð K A basis function w 2 L2(K) defined on a physical element is then constructed from the basis function w the reference element using the standard pull-back transformation:
^ F1 w¼w K
on K:
3.1.2. H1 space In contrast to the space L2, the H1 space requires continuity across element interfaces. The requirement of continuity between elements means that the basis functions must match at vertices and edges. The set of hierarchic basis functions ^ vi are the only naturally splits into sets of functions identified with element vertices, edges and interiors. Vertex functions u non-zero functions at element vertices and therefore must have the same coefficient on all elements containing the vertex.
440
M. Ainsworth, X. Ma / Journal of Computational Physics 231 (2012) 436–453
^ el i are (apart from the vertex functions) the only non-zero functions on element edges and must Similarly, edge functions u ^ Il;m vanish on all shared boundaries have identical coefficients on all elements containing the edge. Finally, interior functions u and their coefficients may be chosen independently on each element. A H1 hierarchic basis functions [17] of order p on the b is defined in Table 1. reference triangle K Similar to the transformation of L2 basis function, the global H1 basis function u defined on a physical element is con^ defined on the reference element using the standard pull-back transformation: structed from the basis function u
^ F1 u¼u on K: K 3.1.3. H(div) space The space H(div) requires the continuity of the normal components across element interfaces. A set of hierarchic basis b is given in Table 2. functions [17] for the vector-valued polynomial space of order p on the reference element K ^ epi varies along the edge ei as a pth order Legendre polynomial. The The normal components of a pth order edge function u ^ Ii;l have vanishing normal components on all three edges (and a interior basis functions can be grouped into two subsets: u Iy^ I^x ^ l;m and u ^ l;m non-zero tangential component on a single edge), and u with vanishing normal and tangential components on all three edges. We refer to the later as the bubble-interior subspace. ^ defined on A global basis function u of H(div) defined on a physical element K is constructed from the basis function u the reference element using the contravariant (Piola) transformation [18]:
u¼
1 ^ Þ F1 ðJ u K detðJK Þ K
on K:
The Piola transformation has the following important property
1 ^ div u detðJK Þ
divu ¼
which is used in the evaluation of the matrices appearing in the system (2). We summarise the distribution and numbers of the degrees of freedom in these spaces in Table 3. 3.2. Assembly of the system Assembly of the global system in the case of non-uniform polynomial approximation using hierarchical basis requires slight modification from the standard finite element sub-assembly. The key observation is that this can be realized easily by the use of appropriate connectivity mappings [19]. Connectivity mappings play a key role in finite element analysis. For instance, they are used to distribute the global solution vector to the elements and in the assembly of the global linear system. The standard finite element sub-assembly procedure for a global matrix M, such as the mass matrix, may be written in the form [19,20]
M¼
X
KK MK K>K ;
ð4Þ
K2P
where MK are element matrices and KK are connectivity mappings. The sub-assembly procedure for a global load vector f may be written in the form
f¼
X
KK f K ;
ð5Þ
K2P
where fK are element load vectors. The distribution of a global solution vector a may be written in the form
aK ¼ Ktop K a;
ð6Þ
where aK are element solution vectors. Table 1 Hierarchic basis functions for H1-conforming finite element space of order p. Degrees of freedom
Basis function
Vertex
^ vi ¼ ^ u ki
Edge 0
i = 0, 1, 2 ^ el 0 ¼ b12 Ll ð^ u n12 Þ
Edge 1 Edge 2
^ el 1 ¼ b20 Ll ð^ u n20 Þ ^ el 2 ¼ b01 Ll ð^ u n01 Þ 06l6p2
Interiors
^ Il;m ¼ b012 Ll ð^ u n20 ÞLm ð^ n01 Þ 0 6 l, m, l + m 6 p 3
M. Ainsworth, X. Ma / Journal of Computational Physics 231 (2012) 436–453
441
Table 2 ^ denotes the gradient with respect to the reference Hierarchic basis functions for H(div)-conforming finite element space, where r element. Degrees of Freedom
Basis Function Lowest order ^^ ^^ ^ e00 ¼ ^ u k1 ð r k2 ðr k2 Þ? ^ k1 Þ ? ? e1 ^ ^ ^ ^ ^ ^ ^ u0 ¼ k2 ðrk0 Þ k0 ðrk2 Þ? ^k ^^ ^1 Þ? k ^0 ð r ^1 ðr ^ e2 ¼ k k0 Þ? u
Edge 0 Edge 1 Edge 2
0
First order ^^ ^^ ^ e10 ¼ ^ k2 Þ? þ ^ k1 Þ ? u k1 ð r k2 ðr ? e1 ^ ^ ^ ^ ^ ^ ^ u1 ¼ k2 ðrk0 Þ þ k0 ðrk2 Þ? ^k ^^ ^1 Þ? þ k ^0 ð r ^1 ðr ^ e2 ¼ k k0 Þ? u
Edge 0 Edge 1 Edge 2
1
(p + 1) th order p 0 ^ ^ e0 ^ ^ e0 ^ epþ1 u ¼ 2pþ1 pþ1 Lp ðn12 Þu1 pþ1 Lp1 ðn12 Þu0
Edge 0 Edge 1
p 1 ^ ^ e1 ^ ^ e1 ^ epþ1 u ¼ 2pþ1 pþ1 Lp ðn20 Þu1 pþ1 Lp1 ðn20 Þu0
Edge 2
p 2 ^ ^ e2 ^ ^ e2 ^ epþ1 u ¼ 2pþ1 pþ1 Lp ðn01 Þu1 pþ1 Lp1 ðn01 Þu0
Interiors
^^ ^ I0;l ¼ b12 ðr u n12 Þ k0 Þ? Ll ð^ ^^ ^ I1;l ¼ b20 ðr u n20 Þ k1 Þ? Ll ð^ ^^ ^ I2;l ¼ b01 ðr u n01 Þ k2 Þ? Ll ð^ l = 0, . . . ,p 2 ^ x ^^ ^ Il;m k1 u ¼ b012 Ll ð^ n20 ÞLm ð^ n01 Þr
I
^ y ^^ ^ l;m k2 u ¼ b012 Ll ð^ n20 ÞLm ð^ n01 Þr
0 6 l, m, l + m 6 p 3
Table 3 Number of degrees of freedom N(c; p) needed on each entity c (i.e. vertex, edge or interior) for the various spaces if the active order of approximation on the entity is p 2 N. (Negative quantities correspond to no basis functions.). Space 1
H =Pp HðdivÞ=Pp Pp L2 =Pp
Vertices
Edges
Interiors
1
p1
(p 2)(p 1)/2
0 0
p+1 0
(p 1)(p + 1) (p + 1)(p + 2)/2
A connectivity mapping is represented by a rectangular matrix of dimension N NK where N is the total number of global degrees of freedom and NK is the number of degrees of freedom on element K. The entries of the matrix are given by
½KK ij ¼
1 if i is the global number of local degree of freedom j; 0
otherwise:
where +1 is chosen if the physical edge and reference edge are oriented in the same sense, otherwise 1 is chosen. In our implementation, the physical edge is oriented from the smaller numbered global vertex to the larger numbered vertex. The reference edge is oriented in anti-clockwise sense as shown in Fig. 1. We emphasize that in order to implement locally non-uniform order approximation, the construction of the connectivity mappings is the only modification of the standard finite element methodology. All other operations such as the evaluation of the element matrices, application of boundary conditions, etc. follow the usual approach. The assignment of global numbers to degrees of freedom in non-uniform approximation is now considered and will be illustrated by considering a simple mesh of four elements
P ¼ fK 0 ; . . . ; K 3 g; shown in Fig. 2. The elements and the desired polynomial orders are given by
K 0 ¼ ½ 0 1 5 ;
p0 ¼ 2;
K 1 ¼ ½ 1 3 5 ;
p1 ¼ 3;
K 2 ¼ ½ 1 2 3 ; K 3 ¼ ½ 3 4 5 ;
p2 ¼ 4; p3 ¼ 5;
where the elements have been oriented in anti-clockwise sense. The second element K1 represents the most common situation whereby an element has neighbouring elements on all edges. The vertices VðK 1 Þ and edges EðK 1 Þ of element K1 are easily identified from the global vertex numbers [1 3 5] of the element:
442
M. Ainsworth, X. Ma / Journal of Computational Physics 231 (2012) 436–453
Fig. 2. Four-element mesh used to illustrate variable order implementation.
Table 4 Vertices and edges of element K1 and local order of approximation determined from using minimum rule. Entity
Vertex labels
Contained within
Local order
v1 v3 v5 e0 e1 e2
[1] [3] [5] [1 3] [3 5] [5 1]
K0, K1, K2, . . . K1, K2, K3, . . . K0, K1, K3, . . . K1, K2 K1, K3 K0, K1
2 3 2 3 3 2
VðK 1 Þ ¼ fv 1 ; v 3 ; v 5 g; EðK 1 Þ ¼ fe0 ; e1 ; e2 g; where the vertices and edges are defined in Table 4. The complete sets V and E of vertices and edges are compiled by forming the unions of such elemental contributions by looping over elements. We need to select a criteria for resolving the active polynomial order on the edges shared by elements of differing nominal polynomial orders. The minimum rule takes the active order of approximation an edge to be the minimum degree of all elements containing the edge [19]. The maximum rule takes the active order of approximation an edge to be the maximum degree of all elements containing the edge. We prefer to use the minimum rule. Table 4 shows the local orders approximation on all vertices and edges of element K1. Later, we shall make use of the maximum rule when performing post-processing and a posteriori error estimation. Once all vertices, edges and interiors have been identified, global numbers are assigned to the degrees of freedom and the connectivity mapping is constructed. This process is performed as follows: 1. For each entity j 2 V [ E [ K, (i.e. each element vertex, edge and interior): (a) Find the local order of approximation pj using the minimum rule. (b) Determine the number of degrees of freedom N(j; pj) needed on the entity by reference to Table 3. (c) Assign global numbers and distribute to all elements containing the entity. 2. Construct the connectivity mapping KK for each element K 2 P: Below, we present in detail the construction of the connectivity mappings for the four element mesh described above in the case of L2, H1 and H(div) conforming approximation. 3.2.1. L2-conforming subspace The basis functions for an L2-conforming subspace given in Eq. (3) consist entirely of interior functions and global numbers are needed only for the interior degrees of freedom, as indicated in Table 3. The local order of approximation on element K1 is p = 3 and Table 3 indicates that 10 degrees of freedom are assigned to element K1. Suppose that the global numbers 116, . . . , 125 are assigned. In theory, the connectivity mapping K1 will be a rectangular matrix of size N 10, where N is the total number of global degrees of freedom with entries defined as in Eq. (4). In practice, this matrix is never explicitly constructed, and instead the information is stored in a more compact form using the representation
RepðK1 Þ ¼ ½116; . . . ; 125: The correspondence between the numbers and the particular basis functions (i.e. the ordering of the functions on the reference element) is a matter of taste, provided that any scheme is used consistently throughout all elements in the mesh.
M. Ainsworth, X. Ma / Journal of Computational Physics 231 (2012) 436–453
443
3.2.2. H1-conforming subspace The basis functions for an H1-conforming subspace given in Table 1 are separated into subsets of basis functions associated with the individual vertices, edges and element interior. The connectivity mapping on element K may be decomposed accordingly, giving
RepðKK Þ ¼ KVK KeK0 KeK1 KeK2 KKK : Table 3 indicates that one degree of freedom per vertex is needed for the space H1 giving a total of three for the element. If the actual numbers assigned are selected to coincide with the vertex numbers, then the compact storage scheme described above gives
Rep KV1 ¼ ½ 1 3 5 Similarly, Table 4 shows that edge e2 has local order of approximation p = 2 and then Table 3 indicates that one global number 9 say, is to be assigned. The local order of element K1 is p = 3 meaning that each edge has two basis functions, while only one freedom number has been assigned. This is to be expected, since the highest order, p = 2 and 3, edge functions must have coefficients equal to zero in order to preserve conformity with lower order elements sharing the edge. This is enforced by assigning a global freedom number of 0 to all such degrees of freedom. The contribution to the connectivity mapping for the element from edge e2 is therefore given by
Rep Ke12 ¼ ½ 9 0 : This process continues until finally it is observed that one interior degree of freedom is needed, 85 say, giving
Rep KK1 ¼ ½85: It is then a simple matter to concatenate the connectivity mappings for each entity to form the element connectivity mapping in the form of a string of integers and zeros. The zeros indicate columns of zeros in the full matrix representation and automatically enforce that the coefficients of the associated basis functions must be zero. This is necessary to ensure conformity since these basis functions on the neighbouring, lower order, elements are not activated. 3.2.3. H(div)-conforming subspace The basis functions for an H(div)-conforming subspace given in Table 2 are separated into subsets of basis functions associated with the edges and element interior and the connectivity mapping takes the form
RepðKK Þ ¼ KeK0 KeK1 KeK2 KKK : The treatment is then identical with that for the H1-case, using the second row of Table 3 to determine the number of degrees of freedom. It will be noticed that in passing from H1 to H(div), the entity with the lowest dimension is dropped at each step. One might wonder what happens in the next step whereby edges would be removed, yielding a mapping of the form
RepðKK Þ ¼ KKK : This simply corresponds to the L2 space considered earlier. Applying the above constructed connectivity mappings, it is a simple matter to assemble the system and distribute the solution by taking the following steps: (1) for all elements K 2 P, calculate local matrices using the nominal polynomial order at the element K. (2) assemble the global matrices and local vector by means of the usual finite element sub-assembly using the Eqs. (4) and (5). (3) after solving the system, distribute the solution to each element using the Eq. (6). Once again, we see the procedure closely follow the standard finite element sub-assembly for the case of the uniform order approximation. The only difference needed for the case of non-uniform order is the determination of the modified connectivity mappings KK. In conclusion, we have shown that non-uniform order approximation can be relatively easily implemented in an existing finite element code. 4. Post-processing procedure It is very common to apply post-processing of mixed FEM to a obtain a smoothed approximation to the potential. Postprocessing schemes for mixed FEM generally take the form of the solution of a local boundary problem for the potential on each element [8,21,22,9]. A post-processing procedure for mixed FEM of non-uniform order has not, to our knowledge, been
444
M. Ainsworth, X. Ma / Journal of Computational Physics 231 (2012) 436–453
considered. Our approach here shares similarities with the above approaches but applies also to the case of non-uniform order. In particular, the post-processing we advocate consists of three steps: (1) Solution of independent local Neumann problems over each element, (2) Inter-element smoothing of the resulting potential, and (3) Solution of independent local Dirichlet problems over each element. 4.1. Solution of local Neumann problem FE Given a mixed finite element approximation uFE on element K of local degree pK, the Neumann post-processing step K ; rK seeks an approximation uK 2 PqK on each element K 2 P such that
ðAgraduK ; gradv ÞK ¼ rFE K ; gradv K R R FE u ¼ K uK ; K K
8v 2 PqK ðKÞ;
ð7Þ
where
qK ¼
2 if pK ¼ 0; pK þ 1 if pK > 0;
The order qK is selected so that the post-processed flux Agrad uK° is a polynomial whose order exceeds that of rFE K by 1. Note that uK° 2 H1(K) and this step occurs only over a single element, which means that the step is inexpensive since it is unnecessary to construct and solve a global matrix. The function uK° can be approximated using the H1 hierarchic basis span{ui : i 2 NK} H1(K) defined in Table 1, by writing uK° in the form
uK ¼
X
ðbK Þi ui :
i2NK
The side constraint on the average value of the potential is applied as follows: Firstly, observe that
Z K
X
uK ¼
ðbK Þi
Z
i2NK
K
ui dx ¼ C>K bK ;
where
Z
ðC K Þi ¼
ui dx:
K
The side constraint then takes the form
C>K bK ¼ lK ; where
lK ¼
Z K
uFE K dx:
The coefficients and lagrange multiplier associated with the constraint are then determined by solving the local system
SK
CK
C>K
0
bK kK
¼
bK
lK
;
where
ðSK Þij ¼ ðbK Þj ¼
R R
K
A$ui $uj dx;
K
rFE K $uj dx:
It is a simple matter to show that this matrix equation is uniquely solvable following the argument given in [16]. Fig. 3(b) shows the functions fuK gK2P obtained in the case of the L-shaped domain problem. 4.2. Inter-element smoothing In general, the piecewise polynomial function uP is discontinuous across element interfaces as indicated in Fig. 3(b). For elements sharing an edge, the same active polynomial order is assigned by the minimum or maximum rules introduced in ~ P 2 H1 ðXÞ on the partition P whose sub-Section 3.2. This helps us to construct a continuous polynomial piecewise function u values at element vertices and edges are given by a weighted averaging of uP as follows:
M. Ainsworth, X. Ma / Journal of Computational Physics 231 (2012) 436–453
445
Fig. 3. Illustration of post-processing procedure. (a) Original mixed FE potential uFE (b) Potential obtained using elementwise Neumann post-processing ~ P (d) H1-conforming potential obtained using elementwise Dirichlet poststep fuK gK2P (c) Smoothed potential using post-processing smoothing step u ~ K gK2P . processing step fu
~P ðxn Þ u
8
K bK ;
j KK j > x K ;
K2P
where jKKj denotes the matrix obtained by replacing all entries in KK by their absolute values. (2) Average the assembled global vector, gives
~ Þ ¼ ðb Þ =xi ; i 2 N: ðb P i P i (3) Distribute the averaged global vector to each element using the standard solution distribution approach (6), gives
~ ¼ KK b ~ : b K P
446
M. Ainsworth, X. Ma / Journal of Computational Physics 231 (2012) 436–453
4.3. Solution of the local Dirichlet problem ~ P as the Dirichlet A Dirichlet problem can be formulated on each element by using the smoothed continuous function u ~ P 2 H1 ðKÞ such that u ~ P jK 2 PqK ðKÞ and condition and seeking u
~P ; gradv ÞK ¼ rFE ðA gradu P ; gradv K ~ P on @K: ~P ¼ u u
8v 2 PqK ðKÞ \ H10 ðKÞ;
To solve this problem, we can re-use the element matrix SK assembled in the initial Neumann post-processing step, after applying the Dirichlet conditions. The above procedure gives an H1-conforming approximation to the potential. The results of mixed FE solution and the three steps of the post-processing procedure are illustrated in Fig. 3 in the case of the problem of finding u 2 H10 ðXÞ such that Du = 1 on the L-shaped domain X :¼ [1, 1]2n[0, 1] [1, 0]. The fact that we deal with only independent problems over each element means that the post-processing is relatively cheap in comparison with the cost of the evaluation of the finite element approximation itself. 5. Application to a posteriori error estimation The following result shows how the post-processed potential defined above can be used to give a computable upper bound on the error in the finite element approximation of the flux: ~ P 2 H1 ðXÞ denote the post-processed potential and let PP : L2 ðXÞ ! M P denote the orthogonal projection, then Theorem 1. Let u 2 FE 2 ~ kr rFE P kA1 6 kA grad uP rP kA1 þ
X
g2K ;
ð8Þ
K2P
where
1
hK
gK ¼ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
p
kmin ðAK Þ
X
kf PP f kK þ
! CðK; cÞkg PN gkc ;
ð9Þ
c2@K\CN
where
CðK; cÞ2 ¼
2hK hK þ Lc ; plc p
ð10Þ
with
Lc ¼ max jx xc j; x2c
and
lc ¼ min jx xc j: x2c
and where xc denotes the vertex of K opposite to the edge c. Moreover, there exists a constant C > 0, independent of any mesh-size, such that
rX ffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 FE ~ P rFE kAgradu g2K : P kA1 6 kr rP kA1 þ C K2P
ð11Þ
Proof. The proof follows closely that of our previous work on lowest order mixed finite element approximation. In particular, we introduce the infinite dimensional problem find (u⁄, r⁄) 2 L2(X) H(div; X) such that n r ¼ n rFE P on CN, and
ðA1 r ; sÞ þ ðu ; divsÞ ¼
R CD
uD n s ds;
ðdivr ; v Þ þ ðPP f ; v Þ ¼ 0 for all (v, s) 2 L2(X) H(div; X). The existence and uniqueness of (u⁄, r⁄) follows the same fashion as the proof of existence and uniqueness for the original problem. FE It will be useful to consider the difference r r⁄ and r rFE P . In particular, both PP f and divrP belong to M P and hence, FE thanks to the definition of uFE P ; rP ,
FE PP f þ divrFE 8v 2 M P P ; v ¼ f þ divrP ; v ¼ 0
447
M. Ainsworth, X. Ma / Journal of Computational Physics 231 (2012) 436–453
⁄ ⁄ FE so that PP f ¼ divrFE P . Likewise, thanks to the definition of (u , r ), we have PP f ¼ divr and therefore, r rP is divergence free. Hence,
r r ; r rFE P
A
1
¼ gradðu u Þ; r rFE ¼ A1 ðr r Þ; r rFE ¼ u u ; div r rFE ¼ 0; P P P
where we have made use of the facts that r = A grad u, r⁄ = A grad u⁄, that u u⁄ = 0 on CD, and that n r ¼ n rFE P on CN. The above identity means that 2 2 FE 2 kr rFE P kA1 ¼ kr r kA1 þ kr rP kA1 ;
and thus, in order to bound the difference r rFE P , it suffices to bound each of the terms appearing on the right hand side. The first term is readily bounded as follows:
kr r k2A1 ¼ ðA1 ðr r Þ; r r Þ ¼ ðgradðu u Þ; r r Þ ¼ ðu u ; divðr r ÞÞ þ ¼ ðu u ; f PP f Þ þ
Z CN
Z
ðu u Þn ðr r Þ
CN
ðu u Þðg n rFE P Þ:
Fig. 4. h adaptive refinement for singular Poisson problem.
4
4
3
3
2
2
1
1
(a) initial order 1
(c) up to order 3
(b) up to order 2
4
4
3
3
2
2
1
1
(d) up to order 4
Fig. 5. Selection of hp adaptively refined meshes for singular Poisson problem. The colour bar indicates the locally polynomial approximation order.
448
M. Ainsworth, X. Ma / Journal of Computational Physics 231 (2012) 436–453 0
10
−1
Error in Energy
10
−2
10
hp "True" error −3
10
hp Error Estimator h1 "True" error h1 Error Estimator
−4
10
1
2
10
10
3
10
4
10
NO. of Degrees of Freedom
Fig. 6. Convergence of h and hp adaptive refinement for singular Poisson problem. 2
1.9
1.8
Effectivity Index
1.7
1.6
1.5
1.4
1.3
1.2
h1 refinement hp refinement
1.1
1 1 10
2
10
3
10
4
10
NO. of Degrees of Freedom
Fig. 7. Effecitvity indices of h and hp adaptive refinement for singular Poisson problem.
Fig. 8. (a) A ten-pole electric motor and (b) Geometry of the ten-pole electric motor.
M. Ainsworth, X. Ma / Journal of Computational Physics 231 (2012) 436–453
Fig. 9. h refined mesh for a ten-pole electric motor.
Observing that,
ðu u ; f PP f ÞK ¼ ðu uK ; f PP f ÞK ; where u = u u⁄ and uK 2 PpK is arbitrary, we obtain from Poincaré’s inequality
ðu u ; f PP f ÞK 6 ku uK kK kf PP f kK 6
hK
p
kgrad ukK kf PP f kK 6
hK
p
1 pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi kgrad ukA;K kf PP f kK : kmin ðAK Þ
449
450
M. Ainsworth, X. Ma / Journal of Computational Physics 231 (2012) 436–453 15 14 13 12 11 10 9 8 7 6 5 4 3 2 1
15 14 13 12 11 10 9 8 7 6 5 4 3 2 1
(a) up to order 2
(b) up to order 3 15 14 13 12 11 10 9 8 7 6 5 4 3 2 1
15 14 13 12 11 10 9 8 7 6 5 4 3 2 1
(c) up to order 6
(d) up to order 9 15 14 13 12 11 10 9 8 7 6 5 4 3 2 1
(e) up to order 12
15 14 13 12 11 10 9 8 7 6 5 4 3 2 1
(f) up to order 15
Fig. 10. hp adaptively refined mesh for a ten-pole electric motor. The color bar indicates the locally polynomial approximation order.
Moreover, for c 2 @K \ CN, using Lemma 11 of [23]
Z c
ðu u Þðg n rFE P Þds ¼
Z c
ðu uK Þðg PN gÞ 6 ku uK kc kg PN gkc 6 CðK; cÞkgrad ukK kg PN gkc
1 6 CðK; cÞ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi kg PN gkc kr r kA1 kmin ðAK Þ and so
kr r k2A1 6
X K2P
or
kgrad ukA1 ;K gK 6
X K2P
!12
g2K
1
kA2 gradðu u Þk ¼
X K2P
!12
g2K
kr r kA1
451
M. Ainsworth, X. Ma / Journal of Computational Physics 231 (2012) 436–453
−3
Absolute error
10
−4
10
hp "True" Error hp Error Estimator h1 "True" Error h1 Error Estimator h2 "True" Error h2 Error Estimator h3 "True" Error h3 Error Estimator h4 "True" Error h4 Error Estimator −5
10
4
5
10
10
NO. of Degrees of Freedom
Fig. 11. Convergence of h and hp adaptive algorithms for a ten-pole electric motor.
2
1.9
1.8
Effectivity index
1.7
1.6
1.5
1.4
1.3
1.2 minimum rule maximum rule
1.1
1 3 10
4
5
10
10
No. of Degrees of Freedom
Fig. 12. Estimated effectivity indices applying minimum and maximum rules by solving a ten-pole electric motor using a hp adaptive algorithm.
kr r k2A1 6
X
g2K ;
K2P
where krkA = (Ar, r) and krkA1 ¼ ðA1 r; rÞ. The bound for the second term is simple:
2 FE FE FE ~ ~ kr rFE P kA1 ¼ r rP ; r A grad uP A1 þ r rP ; A grad uP rP A1 and then note that the first term simplifies to
FE ~ ~ r rFE P ; gradðu uP Þ ¼ div r rP ; u uP þ
Z CD
~ n r rFE P ðu uP Þ;
~ P ¼ 0 on CD. An application of the Cauchy–Schwarz inequality then which vanishes since div r rFE ¼ 0 and since u u P gives 2 FE 2 ~ kr rFE P kA1 6 kAgraduP rP kA1 :
The proof of the lower bound (11) is a simple modification of the arguments given in [16]. As as consequence, we obtain Theorem 1. h
452
M. Ainsworth, X. Ma / Journal of Computational Physics 231 (2012) 436–453
6. Numerical examples We illustrate the post-processing procedure and a posteriori error estimator described in previous section for some representative problems. 6.1. Singular Poisson problem Considering the L-shape domain problem introduced in the last section, we approximate the problem using both h and hp adaptive refinement algorithms with initial mesh shown in Fig. 4(a). The h adaptive refinement algorithm successively refines all elements whose contribution to the error estimator exceeds 30% of the largest local estimator [16]. This leads to the final refined mesh as shown in Fig. 4(b). The hp adaptive refinement algorithm using the ‘‘tagging strategy’’ described in [16] reads as follows: (1) Estimate the error for each element using the post-processing procedure described in Section 3. (2) Flag all the elements whose contribution to the error estimator exceeds 30% of the largest local estimator. (3) Apply h refinement on the flagged elements tagged to the six corner points and p refinement on other flagged elements. (4) Repeat the above steps until the largest local estimator lies within a given tolerance. As shown in Fig. 5, a sequence of hp refined meshes is obtained using the above algorithm. Fig. 6 shows the values of the FE ~ true error in the fluxes kr rFE P kA1 , the estimated error kAgraduP rP kA1 , of the h and hp algorithm. Fig. 7 shows the effectivity indices given by ratio of the estimated error to the true error of the h and hp algorithm. We observe that the estimator indeed provides a guaranteed robust for both h and hp refinement upper bound on the true error in both the h and hp adaptive refinement algorithms. 6.2. Electric motor We illustrate the quality and performance of the estimator for the more challenging example of a ten-pole electric motor shown in Fig. 8. The problem is modelled using the equations of magnetostatics
curl H ¼ J; div B ¼ 0;
ð12Þ
where J is the prescribed current density, where B is the unknown magnetic flux density related to the magnetic field H by the constitutive relation B = lH, and where l is the magnetic permeability of the material. The current density J is unity in the coils and zero value in other domains. The magnetic permeability is given by l = 1 in the coils and the air gap and by l = 5200 in the stator and rotor. Rewrite the second equation in (12) to B = curl u, where u is a scalar magnetic potential, system (12) can be transformed to the equations
H l1 curl u ¼ 0; curl H þ J ¼ 0:
ð13Þ
To reduce this system to the form (1), we set l = A1, J = f, and H = r\. Consequently, the mixed method described in Section 2 can be used to approximate the problem. Considering symmetry, the computations may be carried out on a quarter domain as shown in Fig. 9(a) with homogeneous Dirichlet data on the edges. The geometry of the domain is complex and the material varies significantly in different domains. In particular, we expect the solution to have singularities of differing intensities at each of the re-entrant corners and at material cross-points. It is essential to correctly balance the refinements at each such point if one is to obtain good rates of convergence. The same h adaptive refinement algorithms are applied to the motor problem, using fixed approximation order 1, 2, 3 and 4. Fig. 9(b–e) shows the initial and final h refined meshes. The hp adaptive refinement algorithm applies h refinement on elements tagged to all the corner points and p refinement on other elements. Fig. 10 shows a sequence of hp refined meshes successively refining all elements whose contribution to the error estimator exceeds 30% of the largest local estimator. Then, in Fig. 11, we illustrate the convergence rates of the h and hp adaptive procedures by showing the values of the estimator obtained during the adaptive refinement procedure. Fig. 12 gives the effectivity indices, applying the minimum and maximum rules to approximate the local Neumann problem (7) in the post-processing procedure, by solving the ten-pole electric motor using the hp adaptive algorithm. Once again, we see that the estimator provides very satisfactory guaranteed upper bounds on the true error in both the h and hp adaptive refinement algorithm. The hp algorithm converges fastest compared with the four h adaptive algorithms. 7. Concluding remarks We have described a systematic and complete approach for the mixed finite element approximation of elliptic problems, along with theoretical justification. The approach encompasses a simple post-processing procedure similar to those
M. Ainsworth, X. Ma / Journal of Computational Physics 231 (2012) 436–453
453
proposed for low-order mixed finite elements, and shows that it may be used to produce guaranteed and computable error estimates. For ease of exposition, we have restricted the presentation to the two-dimensional case. However, the same approach is applicable to H(div)-conforming approximation three-dimensions based on existing standard H(div)-conforming elements and existing post-processing schemes in conjunction with the same modification presented here. In particular, Theorem 1 generalizes directly to 3D. References [1] C. Schwab, p- and hp-Finite Element Methods, Numerical Mathematics and Scientific Computation, The Clarendon Press, Oxford University Press, New York, 1998. ISBN 0-19-850390-3. Theory and applications in solid and fluid mechanics. [2] B. Guo, I. Babuka, The h-p version of the finite element method. Part 1 the basic approximation results, Comput. Mech. 1 (1986) 21–41. 10.1007/ BF00298636; . [3] B. Guo, I. Babuka, The h-p version of the finite element method. Part 2 general results and applications, Comput. Mech. 1 (1986) 203–220. 10.1007/ BF00272624; . [4] F. Brezzi, M. Fortin, Mixed and Hybrid Finite Element Methods, Springer-Verlag New York, Inc, NY, USA, 1991. ISBN 0-387-97582-9. [5] M. Ainsworth, K. Pinchedez, hp-approximation theory for BDFM and RT finite elements on quadrilaterals, SIAM J. Numer. Anal. 40 (6) (2002) 2047– 2068 (electronic) (2003). doi:10.1137/S0036142901391128. URL . [6] L. Demkowicz, Computing with hp-adaptive finite elements, Chapman & Hall/CRC Applied Mathematics and Nonlinear Science Series, vol. 1, Chapman & Hall/CRC, Boca Raton, FL, 2007. ISBN 978-1-58488-671-6; 1-58488-671-4. One and two dimensional elliptic and Maxwell problems, With 1 CD-ROM (UNIX). [7] L. Demkowicz, J. Kurtz, D. Pardo, M. Paszyn´ski, W. Rachowicz, A. Zdunek, Computing with hp-adaptive finite elements, Chapman & Hall/CRC Applied Mathematics and Nonlinear Science Series, vol. 2, Chapman & Hall, CRC, Boca Raton, FL, 2008. ISBN 978-1-58488-672-3; 1-58488-672-2. Frontiers: three dimensional elliptic and Maxwell problems with applications. [8] D.N. Arnold, F. Brezzi, Mixed and nonconforming finite element methods: implementation, postprocessing and error estimates, RAIRO Modél Math Anal Numér 19 (1) (1985) 7–32. [9] R. Stenberg, Postprocessing schemes for some mixed finite elements, RAIRO Modél Math Anal Numér 25 (1) (1991) 151–167. [10] A. Alonso, Error estimators for a mixed method, Numer. Math. 74 (4) (1996) 385–395. doi:10.1007/s002110050222. . [11] C. Lovadina, R. Stenberg, Energy norm a posteriori error estimates for mixed finite element methods, Math. Comput. 75 (256) (2006) 1659–1674 (electronic). doi:10.1090/S0025-5718-06-01872-2. URL . [12] D. Braess, R. Verfürth, A posteriori error estimators for the Raviart–Thomas element, SIAM J. Numer. Anal. 33 (6) (1996) 2431–2444. doi:10.1137/ S0036142994264079. URL . [13] C. Carstensen, A posteriori error estimate for the mixed finite element method, Math. Comput. 66 (218) (1997) 465–476. doi:10.1090/S0025-5718-9700837-5. URL . [14] M. Vohralı´k, Unified primal formulation-based a priori and a posteriori error analysis of mixed finite element methods, Math. Comput. 79 (272) (2010) 2001–2032. doi:10.1090/S0025-5718-2010-02375-0. URL . [15] M. Vohralı´k, A posteriori error estimates for lowest-order mixed finite element discretizations of convection–diffusion–reaction equations, SIAM J. Numer. Anal. 45 (4) (2007) 1570–1599 (electronic). doi:10.1137/060653184. URL . [16] M. Ainsworth, A posteriori error estimation for lowest order Raviart–Thomas mixed finite elements, SIAM J. Sci. Comput. 30 (1) (2007/08) 189–204. doi:10.1137/06067331X. URL . [17] M. Ainsworth, J. Coyle, Hierarchic finite element bases on unstructured tetrahedral meshes, Internat. J. Numer. Methods Eng. 58 (14) (2003) 2103– 2130. doi:10.1002/nme.847. URL . [18] J.E. Marsden, T.J.R. Hughes, Mathematical Foundations of Elasticity, Dover Publications Inc., New York, 1994. ISBN 0-486-67865-2. Corrected reprint of the 1983 original. [19] M. Ainsworth, B. Senior, Aspects of an adaptive hp-finite element method: adaptive strategy, conforming approximation and efficient solvers, Comput. Methods Appl. Mech. Eng. 1997; 150(1–4): 65–87. (doi:10.1016/S0045-7825(97)00101-1. Symposium on Advances in Computational Mechanics, vol. 2 (Austin, TX, 1997); URL . [20] G. Kron, A set of principles to interconnect the solutions of physical systems, J. Appl. Phys. 24 (1953) 965–980. [21] L.D. Marini, An inexpensive method for the evaluation of the solution of the lowest order Raviart–Thomas mixed method, SIAM J. Numer. Anal. 22 (3) (1985) 493–496. doi:10.1137/0722029. URL . [22] R. Stenberg, Some new families of finite elements for the Stokes equations, Numer. Math. 56 (8) (1990) 827–838. doi:10.1007/BF01405291. URL . [23] M. Ainsworth, A posteriori error estimation for discontinuous Galerkin finite element approximation, SIAM J. Numer. Anal. 45 (4) (2007) 1777–1798 (electronic). doi:10.1137/060665993. URL .
Contents lists available at SciVerse ScienceDirect
Journal of Computational Physics journal homepage: www.elsevier.com/locate/jcp
Non-uniform order mixed FEM approximation: Implementation, post-processing, computable error bound and adaptivity Mark Ainsworth ⇑,1, Xinhui Ma 2 Mathematics Department, Strathclyde University, 26 Richmond Street, Glasgow G1 1XH, Scotland, United Kingdom
a r t i c l e
i n f o
Article history: Received 29 March 2011 Received in revised form 31 August 2011 Accepted 8 September 2011 Available online 21 September 2011 Keywords: Posteriori error estimation Mixed finite element method Computable error bounds Electromagnetics Flow in porous media Post-processing
a b s t r a c t The present work provides a straightforward and focused set of tools and corresponding theoretical support for the implementation of an adaptive high order finite element code with guaranteed error control for the approximation of elliptic problems in mixed form. The work contains: details of the discretisation using non-uniform order mixed finite elements of arbitrarily high order; a new local post-processing scheme for the primary variable; the use of the post-processing scheme in the derivation of new, fully computable bounds for the error in the flux variable; and, an hp-adaptive refinement strategy based on the a posteriori error estimator. Numerical examples are presented illustrating the results obtained when the procedure is applied to a challenging problem involving a ten-pole electric motor with singularities arising from both geometric features and discontinuities in material properties. The procedure is shown to be capable of producing high accuracy numerical approximations with relatively modest numbers of unknowns. Ó 2011 Elsevier Inc. All rights reserved.
1. Introduction Adaptive high order finite element schemes incorporating both local mesh refinement and variable order approximation have proved to be one of the most effective methods for the numerical treatment of problems governed by elliptic partial differential equations [1–3], and are typically referred to as the hp-version of the finite element method. The solutions of such problems typically exhibit local singularities in the neighbourhood of re-entrant corners, points at which there is a change in type of boundary condition, and at interfaces and cross-points where the material properties undergo a jump change in value. The theory and implementation of hp-FEM is reasonably well-understood and widely practised in the case of H1-conforming finite element approximation. However, many practitioners prefer to use a mixed finite element scheme [4] for the approximation of wide classes of problems arising in fluid flow, electromagnetics and elastostatics where one is often more concerned with the flux or stress than in the potential or displacement. Although hp-FEM has been developed and analysed for such problems and shown to produce exponential rates of convergence for problems with the same kinds of singularities described above [5–7], the practical uptake of such methods is comparatively low in comparison with the interest seen for H1-conforming approximation. In search of reasons for the relatively poor uptake, one may look to the fact that the underlying function spaces are less straightforward to discretise in terms of: the construction of appropriate families of hierarchical basis functions; the assembly of the stiffness matrices
⇑ Corresponding author. 1 2
E-mail address: [email protected] (M. Ainsworth). The work of MA was supported in part by EPSRC under grant ‘‘Numerical Algorithms and Intelligent Software for the Evolving HPC Platform’’. The work of XM was supported by the EPSRC grant (EP/E040993/1) ‘‘Adaptive Numerical Methods for Optoelectronic Devices’’.
0021-9991/$ - see front matter Ó 2011 Elsevier Inc. All rights reserved. doi:10.1016/j.jcp.2011.09.011
M. Ainsworth, X. Ma / Journal of Computational Physics 231 (2012) 436–453
437
and load vectors in the case of non-uniform local orders of approximation; and, the implementation of a fully reliable a posteriori error control and adaptive strategy is more complicated. The present work aims to provide the practitioner with a straightforward and focused set of tools and corresponding theoretical support for handling each of these issues. We begin with formulating a problem in mixed form and describing the details of the discretisation of the problem using non-uniform order mixed finite elements of arbitrarily high order. Full details of the necessary procedure are supplied along with representative examples. We then turn to the issue of post-processing of non-uniform order mixed finite element approximations. Post-processing schemes for low order mixed finite element schemes of uniform order have been investigated by a number of authors [8,9]. Here, we generalise these schemes to the case of high order non-uniform order approximation. The schemes produce a smooth (H1-conforming) approximation to the primal unknown through the successive solution of an elementwise local Neumann problem, followed by a local inter-element smoothing and finally, the solution of an elementwise local Dirichlet problem. The local nature of the post-processing means that the cost is negligible in comparison with that entailed in the computation of the mixed finite element approximation and, in addition, is readily implemented. This post-processing procedure is then exploited in the construction of an a posteriori estimator for the error in the flux variable. Although a number of a posteriori error estimators have been proposed for mixed finite element approximation [10–15] none of these estimators produce computable guaranteed bounds for the error in the case of non-uniform order mixed finite element approximation. We extend the analysis presented in [16] for the case of lowest order mixed finite element approximation to the case of non-uniform order approximation and derive a fully computable bound for the error in the flux variable. The cost of computing the error bound is negligible once the post-processed solution is in hand. Finally, we describe an hp-adaptive refinement strategy based on the a posteriori error estimator. The present paper contains all of the details and theory needed to implement a non-uniform order adaptive mixed finite element scheme for the solution of problems in mixed form along with guaranteed error control. Numerical examples are presented illustrating the results obtained when the procedure is applied to the solution of a simple L-shaped model problem, and for a challenging problem involving a ten-pole electric motor with a number of singularities arising from both geometric features and discontinuities in material properties. The procedure is shown to be capable of producing high accuracy numerical approximations with relatively modest numbers of unknowns. 2. Model problem Consider the following model problem:
r A grad u ¼ 0; div r þ f ¼ 0 u ¼ uD nr¼g
in X ð1Þ
on CD on CN ;
where X is a bounded Lipschitz domain in the plane with outer unit normal n on the polygonal boundary C = CD [ CN split into a relatively open Neumann boundary CN and a closed Dirichlet boundary CD. The data f 2 L2(X), g 2 L2(CN) and uD 2 CðCD Þ are given. We assume uD is piecewise polynomial for simplicity. The space L2(X) consists of functions which are square integrable, while the space H1(X) in addition requires all components of the gradients are square integrable. The matrix-valued function A 2 L1 ðX; R22 Þ is assumed to be symmetric positive definite and, for ease of exposition, is assumed piecewise constant on sub-domains of X. The mixed variational formulation of the problem can be written as: find (u, r) 2 L2(X) H(div; X) such that n r = g on CN and
ðA1 r; sÞ þ ðu; div sÞ ¼
R CD
uD n s ds;
ðdiv r; v Þ þ ðf ; v Þ ¼ 0; R where ðv ; sÞ 2 L2 ðXÞ HD ðdiv; XÞ; ðf ; v Þ ¼ X f v . The space H(div; X) denotes the space of vector fields {s 2 L2(X)2 : div s 2 L2(X)}, which we equip with the weighted inner product defined by
ðr; sÞHðdiv;XÞ ¼ ðr; sÞA1 þ ðdiv r; div sÞ; where ðr; sÞA1 ¼ ðA1 r; sÞ, with the associated norms denoted by k kH(div;X) and k kA1 . When x 2 X we shall use the notation (, )x to denote the integral inner product over x and omit the subscript in the case where x is the physical domain X. The subspace HD(div; X) is defined to be {s 2 H(div; X) : n s = 0 on CN}. For discretization, we consider a family of partitions P of the domain X into the union of non-overlapping, shape-regular triangular elements such that the nonempty intersection of a distinct pair of elements is a single common node or single common edge, and such that CD \ CN is a subset of the nodes in P. Each element is assumed to be shape regular, but the mesh maybe highly refined locally. Each element K is assigned a local polynomial degree pK. The variational problem can be discretized by introducing finite dimensional subspaces M P L2 ðXÞ and Q P Hðdiv; XÞ relative to the partitions P as follows:
438
M. Ainsworth, X. Ma / Journal of Computational Physics 231 (2012) 436–453
MP ¼
v 2 L2 ðXÞ : v jK 2 Pp
K1
8K 2 P ;
and
Q P ¼ s 2 Hðdiv; XÞ : sjK 2 PpK PpK 8K 2 P ; where Pk denotes polynomials of total degree at most k. The space Q P corresponds to the Brezzi–Douglas–Marini finite element [4] in the case pK > 0, whilst the lowest order Raviart–Thomas finite element is used in the case pK = 0. The pair of spaces MP and Q P provides a stable approximation of the mixed problem (1). FE We seek a pair uFE 2 M P Q P satisfying n rFE P ; rP P ¼ PN g on CN
R FE A1 rFE P ; s þ uP ; divs ¼ CD uD n s ds; div rFE P ; v þ ðf ; v Þ ¼ 0 for all ðv ; sÞ 2 M P Q P \ HD ðdiv; XÞ. Here, PNg is the piecewise orthogonal projection onto the space PpK ðcÞ for each edge c CN \ @K and all K 2 P. We introduce hierarchic bases M P ¼ spanfwm : m 2 Mg L2 ðXÞ and Q P \ HD ðdiv; XÞ ¼ spanfun : n 2 Ng HD ðdiv; XÞ, which will be defined in next section. The corresponding approximations to r and u are given respectively by
rFE P ¼ uFE P ¼
P n2N
rn un þ rNP ;
P
m2M
um wm :
where rNP 2 Q P is any function satisfying n rNP ¼ PN g on CN . The coefficients are determined by the solutions of the linear system
P n2N
P
n2N
rn A1 un ; ui þ
P
um ðwm ; divui Þ ¼
m2M
rn ðdivun ; wj Þ ¼
R
X
f wj
R CD
uD n ui A1 rNP ; s 8i 2 N;
8j 2 M;
or in matrix form
Br þ Cu ¼ b >
C r¼f
;
ð2Þ
where, for i, n 2 N and j, m 2 M
Bin ¼ ðA1 un ; ui Þ; C im ¼ ðwm ; divui Þ; R bi ¼ CD uD n ui A1 rNP ; ui ; R fj ¼ X f wj : 3. Basis functions and system assembly in the case of locally non-uniform order approximation In this section, we discuss the details of the calculation needed to assemble the matrices appearing in (2) 3.1. Hierarchic basis The reference domain is taken to be the right triangle
b ¼ fð^x; y ^Þ : 1 < y ^ < 1; 1 < ^x < y ^g K with vertices and edges labeled as shown in Fig. 1. Let ^ k0 ; ^ k1 and ^ k2 denote the usual barycentric, or area coordinates, given by
^k0 ¼ 1 ð^x þ y ^Þ; 2
^k1 ¼ 1 ð^x þ 1Þ; 2
^k2 ¼ 1 ðy ^ þ 1Þ: 2
Each edge is described by a pair of vertex numbers ½ i may be assigned a unique parametrization given by
^nij ¼ ^kj ^ki 2 ½1; 1 ij 2 f01; 12; 20g:
j and oriented from vertex i to j in anti-clockwise sense. An edge
M. Ainsworth, X. Ma / Journal of Computational Physics 231 (2012) 436–453
439
Fig. 1. Notation for reference triangle element and orientation of edges.
For later reference, we introduce bubble functions on each edge and triangle by
bij ¼ ^ki ^kj
on edge ij 2 f01; 12; 20g
and
b012 ¼ ^k0 ^k1 ^k2 : Note that these functions vanish on the boundary of the entity with which they are associated. A hierarchic basis on a reference interval I = [1, 1] can be constructed by a set whose lth member function is a polynomial of degree l = 0, 1, . . .. This set can be used to construct general bases on a triangle and is referred to as the primary basis. One obvious choice for the functions would be the monomials of degree l. However, in practice, the function is taken to be the Legendre polynomial Ll of degree l, defined inductively by the relations
L0 ðsÞ ¼ 1;
L1 ðsÞ ¼ s
and
Llþ1 ðsÞ ¼
2l þ 1 l sLl ðsÞ Ll1 ðsÞ; lþ1 lþ1
l ¼ 1; 2; . . .
The primary basis for an interval I may be used to construct polynomial bases on edges by taking the arguments of the primary basis functions to be edge parametrization nij 2 [1, 1]. The functions defined above are used to build discrete subspaces of L2, H1 and H(div) as follows. 3.1.1. L2 space The space L2 imposes no continuity across element interface. This dramatically simplifies the construction of finite eleb is defined by ^ I of order p on the reference triangle K ment space, and a hierarchic basis functions w l;m
^ I ¼ Ll ð^n20 ÞLm ð^n01 Þ 0 6 l; m; l þ m 6 p: w l;m
ð3Þ
The set of hierarchic basis functions consists entirely of interior functions. The term interior reflects the nature of the degrees of freedom associated with the functions. The absence of inter element continuity requirements for the L2 space means the coefficients of the functions can be chosen independently on all elements, i.e., all degrees of freedom are internal to the element. To transform the local basis on the reference element to the global basis on a physical element, let the mapping from the b to the physical element K be denoted by FK : K b ! K. We assume that FK is a differentiable bijection, and reference element K let JK = DFK denote the Jacobian matrix of the transformation. b Þ defined on ^ 2 L2 ð K A basis function w 2 L2(K) defined on a physical element is then constructed from the basis function w the reference element using the standard pull-back transformation:
^ F1 w¼w K
on K:
3.1.2. H1 space In contrast to the space L2, the H1 space requires continuity across element interfaces. The requirement of continuity between elements means that the basis functions must match at vertices and edges. The set of hierarchic basis functions ^ vi are the only naturally splits into sets of functions identified with element vertices, edges and interiors. Vertex functions u non-zero functions at element vertices and therefore must have the same coefficient on all elements containing the vertex.
440
M. Ainsworth, X. Ma / Journal of Computational Physics 231 (2012) 436–453
^ el i are (apart from the vertex functions) the only non-zero functions on element edges and must Similarly, edge functions u ^ Il;m vanish on all shared boundaries have identical coefficients on all elements containing the edge. Finally, interior functions u and their coefficients may be chosen independently on each element. A H1 hierarchic basis functions [17] of order p on the b is defined in Table 1. reference triangle K Similar to the transformation of L2 basis function, the global H1 basis function u defined on a physical element is con^ defined on the reference element using the standard pull-back transformation: structed from the basis function u
^ F1 u¼u on K: K 3.1.3. H(div) space The space H(div) requires the continuity of the normal components across element interfaces. A set of hierarchic basis b is given in Table 2. functions [17] for the vector-valued polynomial space of order p on the reference element K ^ epi varies along the edge ei as a pth order Legendre polynomial. The The normal components of a pth order edge function u ^ Ii;l have vanishing normal components on all three edges (and a interior basis functions can be grouped into two subsets: u Iy^ I^x ^ l;m and u ^ l;m non-zero tangential component on a single edge), and u with vanishing normal and tangential components on all three edges. We refer to the later as the bubble-interior subspace. ^ defined on A global basis function u of H(div) defined on a physical element K is constructed from the basis function u the reference element using the contravariant (Piola) transformation [18]:
u¼
1 ^ Þ F1 ðJ u K detðJK Þ K
on K:
The Piola transformation has the following important property
1 ^ div u detðJK Þ
divu ¼
which is used in the evaluation of the matrices appearing in the system (2). We summarise the distribution and numbers of the degrees of freedom in these spaces in Table 3. 3.2. Assembly of the system Assembly of the global system in the case of non-uniform polynomial approximation using hierarchical basis requires slight modification from the standard finite element sub-assembly. The key observation is that this can be realized easily by the use of appropriate connectivity mappings [19]. Connectivity mappings play a key role in finite element analysis. For instance, they are used to distribute the global solution vector to the elements and in the assembly of the global linear system. The standard finite element sub-assembly procedure for a global matrix M, such as the mass matrix, may be written in the form [19,20]
M¼
X
KK MK K>K ;
ð4Þ
K2P
where MK are element matrices and KK are connectivity mappings. The sub-assembly procedure for a global load vector f may be written in the form
f¼
X
KK f K ;
ð5Þ
K2P
where fK are element load vectors. The distribution of a global solution vector a may be written in the form
aK ¼ Ktop K a;
ð6Þ
where aK are element solution vectors. Table 1 Hierarchic basis functions for H1-conforming finite element space of order p. Degrees of freedom
Basis function
Vertex
^ vi ¼ ^ u ki
Edge 0
i = 0, 1, 2 ^ el 0 ¼ b12 Ll ð^ u n12 Þ
Edge 1 Edge 2
^ el 1 ¼ b20 Ll ð^ u n20 Þ ^ el 2 ¼ b01 Ll ð^ u n01 Þ 06l6p2
Interiors
^ Il;m ¼ b012 Ll ð^ u n20 ÞLm ð^ n01 Þ 0 6 l, m, l + m 6 p 3
M. Ainsworth, X. Ma / Journal of Computational Physics 231 (2012) 436–453
441
Table 2 ^ denotes the gradient with respect to the reference Hierarchic basis functions for H(div)-conforming finite element space, where r element. Degrees of Freedom
Basis Function Lowest order ^^ ^^ ^ e00 ¼ ^ u k1 ð r k2 ðr k2 Þ? ^ k1 Þ ? ? e1 ^ ^ ^ ^ ^ ^ ^ u0 ¼ k2 ðrk0 Þ k0 ðrk2 Þ? ^k ^^ ^1 Þ? k ^0 ð r ^1 ðr ^ e2 ¼ k k0 Þ? u
Edge 0 Edge 1 Edge 2
0
First order ^^ ^^ ^ e10 ¼ ^ k2 Þ? þ ^ k1 Þ ? u k1 ð r k2 ðr ? e1 ^ ^ ^ ^ ^ ^ ^ u1 ¼ k2 ðrk0 Þ þ k0 ðrk2 Þ? ^k ^^ ^1 Þ? þ k ^0 ð r ^1 ðr ^ e2 ¼ k k0 Þ? u
Edge 0 Edge 1 Edge 2
1
(p + 1) th order p 0 ^ ^ e0 ^ ^ e0 ^ epþ1 u ¼ 2pþ1 pþ1 Lp ðn12 Þu1 pþ1 Lp1 ðn12 Þu0
Edge 0 Edge 1
p 1 ^ ^ e1 ^ ^ e1 ^ epþ1 u ¼ 2pþ1 pþ1 Lp ðn20 Þu1 pþ1 Lp1 ðn20 Þu0
Edge 2
p 2 ^ ^ e2 ^ ^ e2 ^ epþ1 u ¼ 2pþ1 pþ1 Lp ðn01 Þu1 pþ1 Lp1 ðn01 Þu0
Interiors
^^ ^ I0;l ¼ b12 ðr u n12 Þ k0 Þ? Ll ð^ ^^ ^ I1;l ¼ b20 ðr u n20 Þ k1 Þ? Ll ð^ ^^ ^ I2;l ¼ b01 ðr u n01 Þ k2 Þ? Ll ð^ l = 0, . . . ,p 2 ^ x ^^ ^ Il;m k1 u ¼ b012 Ll ð^ n20 ÞLm ð^ n01 Þr
I
^ y ^^ ^ l;m k2 u ¼ b012 Ll ð^ n20 ÞLm ð^ n01 Þr
0 6 l, m, l + m 6 p 3
Table 3 Number of degrees of freedom N(c; p) needed on each entity c (i.e. vertex, edge or interior) for the various spaces if the active order of approximation on the entity is p 2 N. (Negative quantities correspond to no basis functions.). Space 1
H =Pp HðdivÞ=Pp Pp L2 =Pp
Vertices
Edges
Interiors
1
p1
(p 2)(p 1)/2
0 0
p+1 0
(p 1)(p + 1) (p + 1)(p + 2)/2
A connectivity mapping is represented by a rectangular matrix of dimension N NK where N is the total number of global degrees of freedom and NK is the number of degrees of freedom on element K. The entries of the matrix are given by
½KK ij ¼
1 if i is the global number of local degree of freedom j; 0
otherwise:
where +1 is chosen if the physical edge and reference edge are oriented in the same sense, otherwise 1 is chosen. In our implementation, the physical edge is oriented from the smaller numbered global vertex to the larger numbered vertex. The reference edge is oriented in anti-clockwise sense as shown in Fig. 1. We emphasize that in order to implement locally non-uniform order approximation, the construction of the connectivity mappings is the only modification of the standard finite element methodology. All other operations such as the evaluation of the element matrices, application of boundary conditions, etc. follow the usual approach. The assignment of global numbers to degrees of freedom in non-uniform approximation is now considered and will be illustrated by considering a simple mesh of four elements
P ¼ fK 0 ; . . . ; K 3 g; shown in Fig. 2. The elements and the desired polynomial orders are given by
K 0 ¼ ½ 0 1 5 ;
p0 ¼ 2;
K 1 ¼ ½ 1 3 5 ;
p1 ¼ 3;
K 2 ¼ ½ 1 2 3 ; K 3 ¼ ½ 3 4 5 ;
p2 ¼ 4; p3 ¼ 5;
where the elements have been oriented in anti-clockwise sense. The second element K1 represents the most common situation whereby an element has neighbouring elements on all edges. The vertices VðK 1 Þ and edges EðK 1 Þ of element K1 are easily identified from the global vertex numbers [1 3 5] of the element:
442
M. Ainsworth, X. Ma / Journal of Computational Physics 231 (2012) 436–453
Fig. 2. Four-element mesh used to illustrate variable order implementation.
Table 4 Vertices and edges of element K1 and local order of approximation determined from using minimum rule. Entity
Vertex labels
Contained within
Local order
v1 v3 v5 e0 e1 e2
[1] [3] [5] [1 3] [3 5] [5 1]
K0, K1, K2, . . . K1, K2, K3, . . . K0, K1, K3, . . . K1, K2 K1, K3 K0, K1
2 3 2 3 3 2
VðK 1 Þ ¼ fv 1 ; v 3 ; v 5 g; EðK 1 Þ ¼ fe0 ; e1 ; e2 g; where the vertices and edges are defined in Table 4. The complete sets V and E of vertices and edges are compiled by forming the unions of such elemental contributions by looping over elements. We need to select a criteria for resolving the active polynomial order on the edges shared by elements of differing nominal polynomial orders. The minimum rule takes the active order of approximation an edge to be the minimum degree of all elements containing the edge [19]. The maximum rule takes the active order of approximation an edge to be the maximum degree of all elements containing the edge. We prefer to use the minimum rule. Table 4 shows the local orders approximation on all vertices and edges of element K1. Later, we shall make use of the maximum rule when performing post-processing and a posteriori error estimation. Once all vertices, edges and interiors have been identified, global numbers are assigned to the degrees of freedom and the connectivity mapping is constructed. This process is performed as follows: 1. For each entity j 2 V [ E [ K, (i.e. each element vertex, edge and interior): (a) Find the local order of approximation pj using the minimum rule. (b) Determine the number of degrees of freedom N(j; pj) needed on the entity by reference to Table 3. (c) Assign global numbers and distribute to all elements containing the entity. 2. Construct the connectivity mapping KK for each element K 2 P: Below, we present in detail the construction of the connectivity mappings for the four element mesh described above in the case of L2, H1 and H(div) conforming approximation. 3.2.1. L2-conforming subspace The basis functions for an L2-conforming subspace given in Eq. (3) consist entirely of interior functions and global numbers are needed only for the interior degrees of freedom, as indicated in Table 3. The local order of approximation on element K1 is p = 3 and Table 3 indicates that 10 degrees of freedom are assigned to element K1. Suppose that the global numbers 116, . . . , 125 are assigned. In theory, the connectivity mapping K1 will be a rectangular matrix of size N 10, where N is the total number of global degrees of freedom with entries defined as in Eq. (4). In practice, this matrix is never explicitly constructed, and instead the information is stored in a more compact form using the representation
RepðK1 Þ ¼ ½116; . . . ; 125: The correspondence between the numbers and the particular basis functions (i.e. the ordering of the functions on the reference element) is a matter of taste, provided that any scheme is used consistently throughout all elements in the mesh.
M. Ainsworth, X. Ma / Journal of Computational Physics 231 (2012) 436–453
443
3.2.2. H1-conforming subspace The basis functions for an H1-conforming subspace given in Table 1 are separated into subsets of basis functions associated with the individual vertices, edges and element interior. The connectivity mapping on element K may be decomposed accordingly, giving
RepðKK Þ ¼ KVK KeK0 KeK1 KeK2 KKK : Table 3 indicates that one degree of freedom per vertex is needed for the space H1 giving a total of three for the element. If the actual numbers assigned are selected to coincide with the vertex numbers, then the compact storage scheme described above gives
Rep KV1 ¼ ½ 1 3 5 Similarly, Table 4 shows that edge e2 has local order of approximation p = 2 and then Table 3 indicates that one global number 9 say, is to be assigned. The local order of element K1 is p = 3 meaning that each edge has two basis functions, while only one freedom number has been assigned. This is to be expected, since the highest order, p = 2 and 3, edge functions must have coefficients equal to zero in order to preserve conformity with lower order elements sharing the edge. This is enforced by assigning a global freedom number of 0 to all such degrees of freedom. The contribution to the connectivity mapping for the element from edge e2 is therefore given by
Rep Ke12 ¼ ½ 9 0 : This process continues until finally it is observed that one interior degree of freedom is needed, 85 say, giving
Rep KK1 ¼ ½85: It is then a simple matter to concatenate the connectivity mappings for each entity to form the element connectivity mapping in the form of a string of integers and zeros. The zeros indicate columns of zeros in the full matrix representation and automatically enforce that the coefficients of the associated basis functions must be zero. This is necessary to ensure conformity since these basis functions on the neighbouring, lower order, elements are not activated. 3.2.3. H(div)-conforming subspace The basis functions for an H(div)-conforming subspace given in Table 2 are separated into subsets of basis functions associated with the edges and element interior and the connectivity mapping takes the form
RepðKK Þ ¼ KeK0 KeK1 KeK2 KKK : The treatment is then identical with that for the H1-case, using the second row of Table 3 to determine the number of degrees of freedom. It will be noticed that in passing from H1 to H(div), the entity with the lowest dimension is dropped at each step. One might wonder what happens in the next step whereby edges would be removed, yielding a mapping of the form
RepðKK Þ ¼ KKK : This simply corresponds to the L2 space considered earlier. Applying the above constructed connectivity mappings, it is a simple matter to assemble the system and distribute the solution by taking the following steps: (1) for all elements K 2 P, calculate local matrices using the nominal polynomial order at the element K. (2) assemble the global matrices and local vector by means of the usual finite element sub-assembly using the Eqs. (4) and (5). (3) after solving the system, distribute the solution to each element using the Eq. (6). Once again, we see the procedure closely follow the standard finite element sub-assembly for the case of the uniform order approximation. The only difference needed for the case of non-uniform order is the determination of the modified connectivity mappings KK. In conclusion, we have shown that non-uniform order approximation can be relatively easily implemented in an existing finite element code. 4. Post-processing procedure It is very common to apply post-processing of mixed FEM to a obtain a smoothed approximation to the potential. Postprocessing schemes for mixed FEM generally take the form of the solution of a local boundary problem for the potential on each element [8,21,22,9]. A post-processing procedure for mixed FEM of non-uniform order has not, to our knowledge, been
444
M. Ainsworth, X. Ma / Journal of Computational Physics 231 (2012) 436–453
considered. Our approach here shares similarities with the above approaches but applies also to the case of non-uniform order. In particular, the post-processing we advocate consists of three steps: (1) Solution of independent local Neumann problems over each element, (2) Inter-element smoothing of the resulting potential, and (3) Solution of independent local Dirichlet problems over each element. 4.1. Solution of local Neumann problem FE Given a mixed finite element approximation uFE on element K of local degree pK, the Neumann post-processing step K ; rK seeks an approximation uK 2 PqK on each element K 2 P such that
ðAgraduK ; gradv ÞK ¼ rFE K ; gradv K R R FE u ¼ K uK ; K K
8v 2 PqK ðKÞ;
ð7Þ
where
qK ¼
2 if pK ¼ 0; pK þ 1 if pK > 0;
The order qK is selected so that the post-processed flux Agrad uK° is a polynomial whose order exceeds that of rFE K by 1. Note that uK° 2 H1(K) and this step occurs only over a single element, which means that the step is inexpensive since it is unnecessary to construct and solve a global matrix. The function uK° can be approximated using the H1 hierarchic basis span{ui : i 2 NK} H1(K) defined in Table 1, by writing uK° in the form
uK ¼
X
ðbK Þi ui :
i2NK
The side constraint on the average value of the potential is applied as follows: Firstly, observe that
Z K
X
uK ¼
ðbK Þi
Z
i2NK
K
ui dx ¼ C>K bK ;
where
Z
ðC K Þi ¼
ui dx:
K
The side constraint then takes the form
C>K bK ¼ lK ; where
lK ¼
Z K
uFE K dx:
The coefficients and lagrange multiplier associated with the constraint are then determined by solving the local system
SK
CK
C>K
0
bK kK
¼
bK
lK
;
where
ðSK Þij ¼ ðbK Þj ¼
R R
K
A$ui $uj dx;
K
rFE K $uj dx:
It is a simple matter to show that this matrix equation is uniquely solvable following the argument given in [16]. Fig. 3(b) shows the functions fuK gK2P obtained in the case of the L-shaped domain problem. 4.2. Inter-element smoothing In general, the piecewise polynomial function uP is discontinuous across element interfaces as indicated in Fig. 3(b). For elements sharing an edge, the same active polynomial order is assigned by the minimum or maximum rules introduced in ~ P 2 H1 ðXÞ on the partition P whose sub-Section 3.2. This helps us to construct a continuous polynomial piecewise function u values at element vertices and edges are given by a weighted averaging of uP as follows:
M. Ainsworth, X. Ma / Journal of Computational Physics 231 (2012) 436–453
445
Fig. 3. Illustration of post-processing procedure. (a) Original mixed FE potential uFE (b) Potential obtained using elementwise Neumann post-processing ~ P (d) H1-conforming potential obtained using elementwise Dirichlet poststep fuK gK2P (c) Smoothed potential using post-processing smoothing step u ~ K gK2P . processing step fu
~P ðxn Þ u
8
K bK ;
j KK j > x K ;
K2P
where jKKj denotes the matrix obtained by replacing all entries in KK by their absolute values. (2) Average the assembled global vector, gives
~ Þ ¼ ðb Þ =xi ; i 2 N: ðb P i P i (3) Distribute the averaged global vector to each element using the standard solution distribution approach (6), gives
~ ¼ KK b ~ : b K P
446
M. Ainsworth, X. Ma / Journal of Computational Physics 231 (2012) 436–453
4.3. Solution of the local Dirichlet problem ~ P as the Dirichlet A Dirichlet problem can be formulated on each element by using the smoothed continuous function u ~ P 2 H1 ðKÞ such that u ~ P jK 2 PqK ðKÞ and condition and seeking u
~P ; gradv ÞK ¼ rFE ðA gradu P ; gradv K ~ P on @K: ~P ¼ u u
8v 2 PqK ðKÞ \ H10 ðKÞ;
To solve this problem, we can re-use the element matrix SK assembled in the initial Neumann post-processing step, after applying the Dirichlet conditions. The above procedure gives an H1-conforming approximation to the potential. The results of mixed FE solution and the three steps of the post-processing procedure are illustrated in Fig. 3 in the case of the problem of finding u 2 H10 ðXÞ such that Du = 1 on the L-shaped domain X :¼ [1, 1]2n[0, 1] [1, 0]. The fact that we deal with only independent problems over each element means that the post-processing is relatively cheap in comparison with the cost of the evaluation of the finite element approximation itself. 5. Application to a posteriori error estimation The following result shows how the post-processed potential defined above can be used to give a computable upper bound on the error in the finite element approximation of the flux: ~ P 2 H1 ðXÞ denote the post-processed potential and let PP : L2 ðXÞ ! M P denote the orthogonal projection, then Theorem 1. Let u 2 FE 2 ~ kr rFE P kA1 6 kA grad uP rP kA1 þ
X
g2K ;
ð8Þ
K2P
where
1
hK
gK ¼ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
p
kmin ðAK Þ
X
kf PP f kK þ
! CðK; cÞkg PN gkc ;
ð9Þ
c2@K\CN
where
CðK; cÞ2 ¼
2hK hK þ Lc ; plc p
ð10Þ
with
Lc ¼ max jx xc j; x2c
and
lc ¼ min jx xc j: x2c
and where xc denotes the vertex of K opposite to the edge c. Moreover, there exists a constant C > 0, independent of any mesh-size, such that
rX ffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 FE ~ P rFE kAgradu g2K : P kA1 6 kr rP kA1 þ C K2P
ð11Þ
Proof. The proof follows closely that of our previous work on lowest order mixed finite element approximation. In particular, we introduce the infinite dimensional problem find (u⁄, r⁄) 2 L2(X) H(div; X) such that n r ¼ n rFE P on CN, and
ðA1 r ; sÞ þ ðu ; divsÞ ¼
R CD
uD n s ds;
ðdivr ; v Þ þ ðPP f ; v Þ ¼ 0 for all (v, s) 2 L2(X) H(div; X). The existence and uniqueness of (u⁄, r⁄) follows the same fashion as the proof of existence and uniqueness for the original problem. FE It will be useful to consider the difference r r⁄ and r rFE P . In particular, both PP f and divrP belong to M P and hence, FE thanks to the definition of uFE P ; rP ,
FE PP f þ divrFE 8v 2 M P P ; v ¼ f þ divrP ; v ¼ 0
447
M. Ainsworth, X. Ma / Journal of Computational Physics 231 (2012) 436–453
⁄ ⁄ FE so that PP f ¼ divrFE P . Likewise, thanks to the definition of (u , r ), we have PP f ¼ divr and therefore, r rP is divergence free. Hence,
r r ; r rFE P
A
1
¼ gradðu u Þ; r rFE ¼ A1 ðr r Þ; r rFE ¼ u u ; div r rFE ¼ 0; P P P
where we have made use of the facts that r = A grad u, r⁄ = A grad u⁄, that u u⁄ = 0 on CD, and that n r ¼ n rFE P on CN. The above identity means that 2 2 FE 2 kr rFE P kA1 ¼ kr r kA1 þ kr rP kA1 ;
and thus, in order to bound the difference r rFE P , it suffices to bound each of the terms appearing on the right hand side. The first term is readily bounded as follows:
kr r k2A1 ¼ ðA1 ðr r Þ; r r Þ ¼ ðgradðu u Þ; r r Þ ¼ ðu u ; divðr r ÞÞ þ ¼ ðu u ; f PP f Þ þ
Z CN
Z
ðu u Þn ðr r Þ
CN
ðu u Þðg n rFE P Þ:
Fig. 4. h adaptive refinement for singular Poisson problem.
4
4
3
3
2
2
1
1
(a) initial order 1
(c) up to order 3
(b) up to order 2
4
4
3
3
2
2
1
1
(d) up to order 4
Fig. 5. Selection of hp adaptively refined meshes for singular Poisson problem. The colour bar indicates the locally polynomial approximation order.
448
M. Ainsworth, X. Ma / Journal of Computational Physics 231 (2012) 436–453 0
10
−1
Error in Energy
10
−2
10
hp "True" error −3
10
hp Error Estimator h1 "True" error h1 Error Estimator
−4
10
1
2
10
10
3
10
4
10
NO. of Degrees of Freedom
Fig. 6. Convergence of h and hp adaptive refinement for singular Poisson problem. 2
1.9
1.8
Effectivity Index
1.7
1.6
1.5
1.4
1.3
1.2
h1 refinement hp refinement
1.1
1 1 10
2
10
3
10
4
10
NO. of Degrees of Freedom
Fig. 7. Effecitvity indices of h and hp adaptive refinement for singular Poisson problem.
Fig. 8. (a) A ten-pole electric motor and (b) Geometry of the ten-pole electric motor.
M. Ainsworth, X. Ma / Journal of Computational Physics 231 (2012) 436–453
Fig. 9. h refined mesh for a ten-pole electric motor.
Observing that,
ðu u ; f PP f ÞK ¼ ðu uK ; f PP f ÞK ; where u = u u⁄ and uK 2 PpK is arbitrary, we obtain from Poincaré’s inequality
ðu u ; f PP f ÞK 6 ku uK kK kf PP f kK 6
hK
p
kgrad ukK kf PP f kK 6
hK
p
1 pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi kgrad ukA;K kf PP f kK : kmin ðAK Þ
449
450
M. Ainsworth, X. Ma / Journal of Computational Physics 231 (2012) 436–453 15 14 13 12 11 10 9 8 7 6 5 4 3 2 1
15 14 13 12 11 10 9 8 7 6 5 4 3 2 1
(a) up to order 2
(b) up to order 3 15 14 13 12 11 10 9 8 7 6 5 4 3 2 1
15 14 13 12 11 10 9 8 7 6 5 4 3 2 1
(c) up to order 6
(d) up to order 9 15 14 13 12 11 10 9 8 7 6 5 4 3 2 1
(e) up to order 12
15 14 13 12 11 10 9 8 7 6 5 4 3 2 1
(f) up to order 15
Fig. 10. hp adaptively refined mesh for a ten-pole electric motor. The color bar indicates the locally polynomial approximation order.
Moreover, for c 2 @K \ CN, using Lemma 11 of [23]
Z c
ðu u Þðg n rFE P Þds ¼
Z c
ðu uK Þðg PN gÞ 6 ku uK kc kg PN gkc 6 CðK; cÞkgrad ukK kg PN gkc
1 6 CðK; cÞ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi kg PN gkc kr r kA1 kmin ðAK Þ and so
kr r k2A1 6
X K2P
or
kgrad ukA1 ;K gK 6
X K2P
!12
g2K
1
kA2 gradðu u Þk ¼
X K2P
!12
g2K
kr r kA1
451
M. Ainsworth, X. Ma / Journal of Computational Physics 231 (2012) 436–453
−3
Absolute error
10
−4
10
hp "True" Error hp Error Estimator h1 "True" Error h1 Error Estimator h2 "True" Error h2 Error Estimator h3 "True" Error h3 Error Estimator h4 "True" Error h4 Error Estimator −5
10
4
5
10
10
NO. of Degrees of Freedom
Fig. 11. Convergence of h and hp adaptive algorithms for a ten-pole electric motor.
2
1.9
1.8
Effectivity index
1.7
1.6
1.5
1.4
1.3
1.2 minimum rule maximum rule
1.1
1 3 10
4
5
10
10
No. of Degrees of Freedom
Fig. 12. Estimated effectivity indices applying minimum and maximum rules by solving a ten-pole electric motor using a hp adaptive algorithm.
kr r k2A1 6
X
g2K ;
K2P
where krkA = (Ar, r) and krkA1 ¼ ðA1 r; rÞ. The bound for the second term is simple:
2 FE FE FE ~ ~ kr rFE P kA1 ¼ r rP ; r A grad uP A1 þ r rP ; A grad uP rP A1 and then note that the first term simplifies to
FE ~ ~ r rFE P ; gradðu uP Þ ¼ div r rP ; u uP þ
Z CD
~ n r rFE P ðu uP Þ;
~ P ¼ 0 on CD. An application of the Cauchy–Schwarz inequality then which vanishes since div r rFE ¼ 0 and since u u P gives 2 FE 2 ~ kr rFE P kA1 6 kAgraduP rP kA1 :
The proof of the lower bound (11) is a simple modification of the arguments given in [16]. As as consequence, we obtain Theorem 1. h
452
M. Ainsworth, X. Ma / Journal of Computational Physics 231 (2012) 436–453
6. Numerical examples We illustrate the post-processing procedure and a posteriori error estimator described in previous section for some representative problems. 6.1. Singular Poisson problem Considering the L-shape domain problem introduced in the last section, we approximate the problem using both h and hp adaptive refinement algorithms with initial mesh shown in Fig. 4(a). The h adaptive refinement algorithm successively refines all elements whose contribution to the error estimator exceeds 30% of the largest local estimator [16]. This leads to the final refined mesh as shown in Fig. 4(b). The hp adaptive refinement algorithm using the ‘‘tagging strategy’’ described in [16] reads as follows: (1) Estimate the error for each element using the post-processing procedure described in Section 3. (2) Flag all the elements whose contribution to the error estimator exceeds 30% of the largest local estimator. (3) Apply h refinement on the flagged elements tagged to the six corner points and p refinement on other flagged elements. (4) Repeat the above steps until the largest local estimator lies within a given tolerance. As shown in Fig. 5, a sequence of hp refined meshes is obtained using the above algorithm. Fig. 6 shows the values of the FE ~ true error in the fluxes kr rFE P kA1 , the estimated error kAgraduP rP kA1 , of the h and hp algorithm. Fig. 7 shows the effectivity indices given by ratio of the estimated error to the true error of the h and hp algorithm. We observe that the estimator indeed provides a guaranteed robust for both h and hp refinement upper bound on the true error in both the h and hp adaptive refinement algorithms. 6.2. Electric motor We illustrate the quality and performance of the estimator for the more challenging example of a ten-pole electric motor shown in Fig. 8. The problem is modelled using the equations of magnetostatics
curl H ¼ J; div B ¼ 0;
ð12Þ
where J is the prescribed current density, where B is the unknown magnetic flux density related to the magnetic field H by the constitutive relation B = lH, and where l is the magnetic permeability of the material. The current density J is unity in the coils and zero value in other domains. The magnetic permeability is given by l = 1 in the coils and the air gap and by l = 5200 in the stator and rotor. Rewrite the second equation in (12) to B = curl u, where u is a scalar magnetic potential, system (12) can be transformed to the equations
H l1 curl u ¼ 0; curl H þ J ¼ 0:
ð13Þ
To reduce this system to the form (1), we set l = A1, J = f, and H = r\. Consequently, the mixed method described in Section 2 can be used to approximate the problem. Considering symmetry, the computations may be carried out on a quarter domain as shown in Fig. 9(a) with homogeneous Dirichlet data on the edges. The geometry of the domain is complex and the material varies significantly in different domains. In particular, we expect the solution to have singularities of differing intensities at each of the re-entrant corners and at material cross-points. It is essential to correctly balance the refinements at each such point if one is to obtain good rates of convergence. The same h adaptive refinement algorithms are applied to the motor problem, using fixed approximation order 1, 2, 3 and 4. Fig. 9(b–e) shows the initial and final h refined meshes. The hp adaptive refinement algorithm applies h refinement on elements tagged to all the corner points and p refinement on other elements. Fig. 10 shows a sequence of hp refined meshes successively refining all elements whose contribution to the error estimator exceeds 30% of the largest local estimator. Then, in Fig. 11, we illustrate the convergence rates of the h and hp adaptive procedures by showing the values of the estimator obtained during the adaptive refinement procedure. Fig. 12 gives the effectivity indices, applying the minimum and maximum rules to approximate the local Neumann problem (7) in the post-processing procedure, by solving the ten-pole electric motor using the hp adaptive algorithm. Once again, we see that the estimator provides very satisfactory guaranteed upper bounds on the true error in both the h and hp adaptive refinement algorithm. The hp algorithm converges fastest compared with the four h adaptive algorithms. 7. Concluding remarks We have described a systematic and complete approach for the mixed finite element approximation of elliptic problems, along with theoretical justification. The approach encompasses a simple post-processing procedure similar to those
M. Ainsworth, X. Ma / Journal of Computational Physics 231 (2012) 436–453
453
proposed for low-order mixed finite elements, and shows that it may be used to produce guaranteed and computable error estimates. For ease of exposition, we have restricted the presentation to the two-dimensional case. However, the same approach is applicable to H(div)-conforming approximation three-dimensions based on existing standard H(div)-conforming elements and existing post-processing schemes in conjunction with the same modification presented here. In particular, Theorem 1 generalizes directly to 3D. References [1] C. Schwab, p- and hp-Finite Element Methods, Numerical Mathematics and Scientific Computation, The Clarendon Press, Oxford University Press, New York, 1998. ISBN 0-19-850390-3. Theory and applications in solid and fluid mechanics. [2] B. Guo, I. Babuka, The h-p version of the finite element method. Part 1 the basic approximation results, Comput. Mech. 1 (1986) 21–41. 10.1007/ BF00298636; . [3] B. Guo, I. Babuka, The h-p version of the finite element method. Part 2 general results and applications, Comput. Mech. 1 (1986) 203–220. 10.1007/ BF00272624; . [4] F. Brezzi, M. Fortin, Mixed and Hybrid Finite Element Methods, Springer-Verlag New York, Inc, NY, USA, 1991. ISBN 0-387-97582-9. [5] M. Ainsworth, K. Pinchedez, hp-approximation theory for BDFM and RT finite elements on quadrilaterals, SIAM J. Numer. Anal. 40 (6) (2002) 2047– 2068 (electronic) (2003). doi:10.1137/S0036142901391128. URL . [6] L. Demkowicz, Computing with hp-adaptive finite elements, Chapman & Hall/CRC Applied Mathematics and Nonlinear Science Series, vol. 1, Chapman & Hall/CRC, Boca Raton, FL, 2007. ISBN 978-1-58488-671-6; 1-58488-671-4. One and two dimensional elliptic and Maxwell problems, With 1 CD-ROM (UNIX). [7] L. Demkowicz, J. Kurtz, D. Pardo, M. Paszyn´ski, W. Rachowicz, A. Zdunek, Computing with hp-adaptive finite elements, Chapman & Hall/CRC Applied Mathematics and Nonlinear Science Series, vol. 2, Chapman & Hall, CRC, Boca Raton, FL, 2008. ISBN 978-1-58488-672-3; 1-58488-672-2. Frontiers: three dimensional elliptic and Maxwell problems with applications. [8] D.N. Arnold, F. Brezzi, Mixed and nonconforming finite element methods: implementation, postprocessing and error estimates, RAIRO Modél Math Anal Numér 19 (1) (1985) 7–32. [9] R. Stenberg, Postprocessing schemes for some mixed finite elements, RAIRO Modél Math Anal Numér 25 (1) (1991) 151–167. [10] A. Alonso, Error estimators for a mixed method, Numer. Math. 74 (4) (1996) 385–395. doi:10.1007/s002110050222. . [11] C. Lovadina, R. Stenberg, Energy norm a posteriori error estimates for mixed finite element methods, Math. Comput. 75 (256) (2006) 1659–1674 (electronic). doi:10.1090/S0025-5718-06-01872-2. URL . [12] D. Braess, R. Verfürth, A posteriori error estimators for the Raviart–Thomas element, SIAM J. Numer. Anal. 33 (6) (1996) 2431–2444. doi:10.1137/ S0036142994264079. URL . [13] C. Carstensen, A posteriori error estimate for the mixed finite element method, Math. Comput. 66 (218) (1997) 465–476. doi:10.1090/S0025-5718-9700837-5. URL . [14] M. Vohralı´k, Unified primal formulation-based a priori and a posteriori error analysis of mixed finite element methods, Math. Comput. 79 (272) (2010) 2001–2032. doi:10.1090/S0025-5718-2010-02375-0. URL . [15] M. Vohralı´k, A posteriori error estimates for lowest-order mixed finite element discretizations of convection–diffusion–reaction equations, SIAM J. Numer. Anal. 45 (4) (2007) 1570–1599 (electronic). doi:10.1137/060653184. URL . [16] M. Ainsworth, A posteriori error estimation for lowest order Raviart–Thomas mixed finite elements, SIAM J. Sci. Comput. 30 (1) (2007/08) 189–204. doi:10.1137/06067331X. URL . [17] M. Ainsworth, J. Coyle, Hierarchic finite element bases on unstructured tetrahedral meshes, Internat. J. Numer. Methods Eng. 58 (14) (2003) 2103– 2130. doi:10.1002/nme.847. URL . [18] J.E. Marsden, T.J.R. Hughes, Mathematical Foundations of Elasticity, Dover Publications Inc., New York, 1994. ISBN 0-486-67865-2. Corrected reprint of the 1983 original. [19] M. Ainsworth, B. Senior, Aspects of an adaptive hp-finite element method: adaptive strategy, conforming approximation and efficient solvers, Comput. Methods Appl. Mech. Eng. 1997; 150(1–4): 65–87. (doi:10.1016/S0045-7825(97)00101-1. Symposium on Advances in Computational Mechanics, vol. 2 (Austin, TX, 1997); URL . [20] G. Kron, A set of principles to interconnect the solutions of physical systems, J. Appl. Phys. 24 (1953) 965–980. [21] L.D. Marini, An inexpensive method for the evaluation of the solution of the lowest order Raviart–Thomas mixed method, SIAM J. Numer. Anal. 22 (3) (1985) 493–496. doi:10.1137/0722029. URL . [22] R. Stenberg, Some new families of finite elements for the Stokes equations, Numer. Math. 56 (8) (1990) 827–838. doi:10.1007/BF01405291. URL . [23] M. Ainsworth, A posteriori error estimation for discontinuous Galerkin finite element approximation, SIAM J. Numer. Anal. 45 (4) (2007) 1777–1798 (electronic). doi:10.1137/060665993. URL .
