A new discretization for mth-Laplace equations with arbitrary polynomial degrees
aa r X i v : . [ m a t h . NA ] J u l A new discretization for m th-Laplace equationswith arbitrary polynomial degrees ∗ M. Schedensack † Abstract
This paper introduces new mixed formulations and discretizations for m th-Laplaceequations of the form ( − m ∆ m u = f for arbitrary m = 1 , , , . . . based on novelHelmholtz-type decompositions for tensor-valued functions. The new discretizationsallow for ansatz spaces of arbitrary polynomial degree and the lowest-order choicecoincides with the non-conforming FEMs of Crouzeix and Raviart for m = 1 and ofMorley for m = 2 . Since the derivatives are directly approximated, the lowest-orderdiscretizations consist of piecewise affine and piecewise constant functions for any m = 1 , , . . . Moreover, a uniform implementation for arbitrary m is possible. Besidesthe a priori and a posteriori analysis, this paper proves optimal convergence rates foradaptive algorithms for the new discretizations. Keywords m th-Laplace equation, polyharmonic equation, non-conforming FEM, mixed FEM,adaptive FEM, optimality AMS subject classification
This paper considers m th-Laplace equations of the form ( − m ∆ m u = f (1.1)for arbitrary m = 1 , , , . . . Standard conforming FEMs require ansatz spaces in H m (Ω) .To circumvent those high regularity requirements and resulting complicated finite elements,non-standard methods are of high interest [Mor68, EGH +
02, Bre12, GN11]. The novelHelmholtz decomposition of this paper decomposes any (tensor-valued) L function in an m th derivative and a symmetric part of a Curl. Given a tensor-valued function ϕ whichsatisfies − div m ϕ = f in the weak sense, the L projection of ϕ to the space D m H m (Ω) of m th derivatives then coincides with the m th derivative of the exact solution of (1.1) (seeTheorem 5.1 below). This results in novel mixed formulations and discretizations for (1.1).This approach generalises the discretizations of [Sch15, Sch16] from m = 1 to m ≥ .The direct approximation of D m u instead of u enables low order discretizations; onlyfirst derivatives appear in the symmetric part of the Curl and so the lowest order approachonly requires piecewise affine functions for any m . In contrast to that, even interior penaltymethods require piecewise quadratic [Bre12] resp. piecewise cubic [GN11] functions for m = 2 resp. m = 3 . Mnemonic diagrams in Figure 1 illustrate lowest-order standardconforming FEMs from [Žen70] and the lowest-order novel FEMs proposed in this work for m = 2 , . Since the proposed new FEMs differ only in the number of components in the ∗ This work was supported by the Berlin Mathematical School. † Institut für Numerische Simulation, Universität Bonn, Wegelerstraße 6, D-53115 Bonn, Germany = 2 3 × × m = 3 4 × × Figure 1: Lowest order standard conforming [Cia78, Žen70] and novel FEMs for the problem ( − m ∆ m u = f for m = 2 , . ansatz spaces, an implementation of one single program, which runs for arbitrary order, ispossible. In particular, the system matrices are obtained by integration of standard FEMbasis functions.For m = 1 , and the lowest polynomial degree in the ansatz spaces, discrete Helmholtzdecompositions of [AF89, CGH14] prove that the discrete solutions are piecewise gradients(resp. Hessians) of Crouzeix-Raviart [CR73] (resp. Morley [Mor68]) finite element func-tions and therefore the new discretizations can be regarded as a generalization of thosenon-conforming FEMs to higher polynomial degrees and higher-order problems. The gen-eralization of [WX13] of the non-conforming Crouzeix-Raviart and Morley FEMs to m ≥ is restricted to a space dimension ≥ m .In the context of the novel (mixed) formulations, the discretizations appear to be con-forming. The new generalization to higher polynomial degrees proposed in this paperappears to be natural in the sense that the inherent properties of the lowest order dis-cretization carry over to higher polynomial ansatz spaces, namely an inf-sup condition, theconformity of the method, and a crucial projection property (also known as integral meanproperty of the non-conforming interpolation operator).Besides the a priori and a posteriori error analysis, this paper proves optimal conver-gence rates for an adaptive algorithm, which are also observed in the numerical experimentsfrom Section 7.The remaining parts of this paper are organised as follows. Section 2 introduces somenotation while some preliminary results are proved in Section 3. The proposed discretiza-tion of (1.1) in Section 5 is based on a novel Helmholtz decomposition for higher derivativeswhich is stated and proved in Section 4. Section 6 introduces an adaptive algorithm andproves optimal convergence rates. Section 7 concludes the paper with numerical experi-ments on fourth- and sixth-order problems.Throughout this paper, let Ω ⊆ R be a bounded, polygonal, simply connected Lip-schitz domain. Standard notation on Lebesgue and Sobolev spaces and their norms isemployed with L scalar product ( • , • ) L (Ω) . Given a Hilbert space X , let L (Ω; X ) resp. H k (Ω; X ) denote the space of functions with values in X whose components are in L (Ω) resp. H k (Ω) . The space of infinitely differentiable functions reads C ∞ (Ω) and the subspaceof functions with compact support in Ω is denoted with C ∞ c (Ω) . The piecewise action ofdifferential operators is denoted with a subscript NC . The formula A . B represents aninequality A ≤ CB for some mesh-size independent, positive generic constant C ; A ≈ B abbreviates A . B . A . By convention, all generic constants C ≈ do neither dependon the mesh-size nor on the level of a triangulation but may depend on the fixed coarsetriangulation T and its interior angles. 2 Notation
This section introduces notation related to higher-order tensors and tensor-valued functionsand triangulations.Define the set of ℓ -tensors over R by X ( ℓ ) := ( R for ℓ = 0 , Q ℓj =1 R = R × · · · × R ∼ = R ℓ for ℓ ≥ and let S ℓ := { σ : { , . . . , ℓ } → { , . . . , ℓ } | σ is bijective } denote the symmetric group, i.e.,the set of all permutations of (1 , . . . , ℓ ) . Define the set of symmetric tensors S ( ℓ ) ⊆ X ( ℓ ) by S ( ℓ ) := { A ∈ X ( ℓ ) | ∀ ( j , . . . , j ℓ ) ∈ { , } ℓ ∀ σ ∈ S ℓ : A j ,...,j ℓ = A j σ (1) ,...,j σ ( ℓ ) } . The symmetric part sym A ∈ S ( ℓ ) of a tensor A ∈ X ( ℓ ) is defined by (sym A ) j ,...,j ℓ := (card( S ℓ )) − X σ ∈ S ℓ A j σ (1) ,...,j σ ( ℓ ) for all ( j , . . . , j ℓ ) ∈ { , } ℓ , where card( M ) denotes the number of elements in a set M .For ℓ = 2 , the set S (2) coincides with the set of symmetric × matrices, while for ℓ = 3 ,the tensors A ∈ S (3) consist of the four different components A , A = A = A , A = A = A , and A . Given ℓ -tensors A, B ∈ X ( ℓ ) and a vector q ∈ R , definethe scalar product A : B ∈ R and the dot product A · q ∈ X ( ℓ − by A : B := X ( j ,...,j ℓ ) ∈{ , } ℓ A j ,...,j ℓ B j ,...,j ℓ , ( A · q ) j ,...,j ℓ − := A j ,...,j ℓ − , q + A j ,...,j ℓ − , q for all ( j , . . . , j ℓ − ) ∈ { , } ℓ − . The following definition summarizes some differentialoperators. Recall that, for a Hilbert space X , the space H (Ω; X ) (resp. L (Ω; X ) ) denotesthe space of H (resp. L ) functions with components in X . Definition 1 (differential operators) . Let v ∈ H ℓ (Ω) and σ ∈ H (Ω; X ( ℓ )) and define p : { , } → { , } by p (1) = 2 and p (2) = 1 . Define the ℓ th derivative D ℓ v ∈ L (Ω; X ( ℓ )) of v , the derivative Dσ ∈ L (Ω; X ( ℓ + 1)) , the divergence div σ ∈ L (Ω; X ( ℓ − , the Curl, Curl σ ∈ L (Ω; X ( ℓ + 1)) , and the curl, curl σ ∈ L (Ω; X ( ℓ − by ( D ℓ v ) j ,...,j ℓ := ∂ ℓ v/ ( ∂x j . . . ∂x j ℓ ) , ( Dσ ) j ,...,j ℓ +1 := ∂σ j ,...,j ℓ /∂x j ℓ +1 , (Curl σ ) j ,...,j ℓ +1 := ( − j ℓ +1 ∂σ j ,...,j ℓ /∂x p ( j ℓ +1 ) , (div σ ) j ,...,j ℓ − := ∂σ j ,...,j ℓ − , /∂x + ∂σ j ,...,j ℓ − , /∂x , (curl σ ) j ,...,j ℓ − := − ∂σ j ,...,j ℓ − , /∂x + ∂σ j ,...,j ℓ − , /∂x for ( j , . . . , j ℓ +1 ) ∈ { , } ℓ +1 .For ℓ = 2 , these definitions coincide with the row-wise application of D , div , Curl , and curl . The L scalar product ( • , • ) L (Ω) of tensor-valued functions f, g : Ω → X ( ℓ ) is definedby ( f, g ) L (Ω) := ´ Ω f : g dx . Given ψ ∈ L (Ω; X ( ℓ )) such that there exists g ∈ L (Ω) with ( ψ, D ℓ v ) L (Ω) = ( − ℓ ( g, v ) L (Ω) for all v ∈ H ℓ (Ω) , ℓ th order divergence div ℓ ψ := g of ψ . The space H (div ℓ , Ω) ⊆ L (Ω; X ( ℓ )) isdefined by H (div ℓ , Ω) := { ψ ∈ L (Ω; X ( ℓ )) | div ℓ ψ ∈ L (Ω) } . Define furthermore for k ≥ ℓH (div ℓ , Ω; X ( k )) := (cid:26) ψ ∈ L (Ω; X ( k )) (cid:12)(cid:12)(cid:12)(cid:12) ∀ ( j , . . . , j k − ℓ ) ∈ { , } k − ℓ : ψ j ,...,j k − ℓ , • ∈ H (div ℓ , Ω) (cid:27) . Remark . Note that the existence of the ℓ th weak divergence does not imply the exis-tence of any k -th divergence for ≤ k ≤ ℓ , e.g., H (div , Ω; X ( ℓ )) H (div ℓ , Ω) for ℓ > .A shape-regular triangulation T of a bounded, polygonal, open Lipschitz domain Ω ⊆ R is a set of closed triangles T ∈ T such that Ω = S T and any two distinct trianglesare either disjoint or share exactly one common edge or one vertex. Let E ( T ) denote theedges of a triangle T and E := E ( T ) := S T ∈ T E ( T ) the set of edges in T . Any edge E ∈ E is associated with a fixed orientation of the unit normal ν E on E (and τ E = (0 , −
1; 1 , ν E denotes the unit tangent on E ). On the boundary, ν E is the outer unit normal of Ω ,while for interior edges E ∂ Ω , the orientation is fixed through the choice of the triangles T + ∈ T and T − ∈ T with E = T + ∩ T − and ν E := ν T + | E is the outer normal of T + on E .In this situation, [ v ] E := v | T + − v | T − denotes the jump across E . For an edge E ⊆ ∂ Ω onthe boundary, the jump across E reads [ v ] E := v . For T ∈ T and X ⊆ X ( ℓ ) , let P k ( T ; X ) := { v : T → X | each component of v is a polynomial of total degree ≤ k } ; P k ( T ; X ) := { v : Ω → X | ∀ T ∈ T : v | T ∈ P k ( T ; X ) } denote the set of piecewise polynomials and P k ( T ) := P k ( T ; R ) . Given a subspace X ⊆ L (Ω; X ( ℓ )) , let Π X : L (Ω; X ( ℓ )) → X denote the L projection onto X and let Π k abbreviate Π P k ( T ; X ( ℓ )) . Given a triangle T ∈ T , let h T := (meas ( T )) / denote the squareroot of the area of T and let h T ∈ P ( T ) denote the piecewise constant mesh-size with h T | T := h T for all T ∈ T . For a set of triangles M ⊆ T , let k • k M abbreviate k • k M := s X T ∈ M k • k L ( T ) . The main result of this section is Theorem 3.2, which proves that k sym Curl •k L (Ω) definesa norm on the space Y defined in (3.5) below and can, thus, be viewed as a generalizedKorn inequality. The following theorem is used in the proof of Theorem 3.2. Recall thedefinition of the Curl and the symmetric part of a tensor from Section 2. Theorem 3.1.
Any γ ∈ H (Ω; S ( m − satisfies k Curl γ k L (Ω) . k sym Curl γ k L (Ω) + k γ k L (Ω) . Proof.
The proof is subdivided in three steps.
Step 1.
Let ≤ k ≤ m and j ( k ) = ( j , . . . , j m ) ∈ { , } m with j ℓ = 1 for all ℓ ∈ { , . . . , k } and j ℓ = 2 for all ℓ ∈ { k + 1 , . . . , m } , i.e., j ( k ) = (1 , . . . , | {z } k , , . . . , | {z } m − k ) . sym and Curl reads (sym Curl γ ) j ( k ) = card( S m ) − X σ ∈ S m ( − j σ ( m ) ∂∂x p ( j σ ( m ) ) γ j σ (1) ,...,j σ ( m − . (3.1)Let j ( k ) := ( j , . . . , j m − ) ∈ { , } m − be the multi-index with the same number of onesand the number of twos reduced by one and j ( k ) := ( j , . . . , j m ) ∈ { , } m − the multi-index with the same number of twos and the number of ones reduced by one, i.e., j ( k ) = (1 , . . . , | {z } k , , . . . , | {z } m − k − ) and j ( k ) = (1 , . . . , | {z } k − , , . . . , | {z } m − k ) . The symmetry of γ implies that γ j σ (1) ,...,j σ ( m − = γ j ( k ) if j σ ( m ) = 1 and γ j σ (1) ,...,j σ ( m − = γ j ( k ) if j σ ( m ) = 2 . Since the number of permutations σ ∈ S m such that j σ ( m ) = 1 is k card( S m − ) and the number of permutations σ ∈ S m such that j σ ( m ) = 2 is ( m − k ) card( S m − ) and since card( S m ) = m ! and card( S m − ) = ( m − , this implies that(3.1) equals (sym Curl γ ) j ( k ) = card( S m − )card( S m ) ( m − k ) ∂γ j ( k ) ∂x − k ∂γ j ( k ) ∂y ! = m − km ∂γ j ( k ) ∂x − km ∂γ j ( k ) ∂y . (3.2) Step 2.
This step applies [Neč67, Chap. 3, Thm. 7.6] and [Neč67, Chap. 3, Thm. 7.8] tooperators N k defined below. Step 3 then proves a relation between these operators andthe operator sym Curl .Define for k ∈ { , . . . , m + 1 } , s ∈ { , . . . , m } , and a multi-index κ ∈ { (1 , , (0 , } a k,s,κ := − ( k − /m if s = k − and κ = (0 , , ( m − k + 1) /m if s = k and κ = (1 , , else . Furthermore, define for ξ ∈ R N ks ξ := X κ =(1 , , (0 , a k,s,κ ξ κ = − ( k − ξ /m if s = k − , ( m − k + 1) ξ /m if s = k, elsewith the multi-index notation ξ κ = ξ κ ξ κ . Then the matrix ( N ks ξ ) ≤ k ≤ ( m +1)1 ≤ s ≤ m reads m mξ . . . − ξ ( m − ξ . . . − ξ ( m − ξ . . . ... ...... . . . . . . ...... ... . . . − ( m − ξ ξ . . . − ( m − ξ ξ . . . − mξ ∈ R ( m +1) × m . ξ = 0 , the columns of this matrix are linear independent. Define the operators ( N k ) k =1 ,...,m +1 , N k : H (Ω; R m ) → L (Ω) , by N k v := m X s =1 X κ =(1 , , (0 , a k,s,κ D κ v s . Then, the combination of [Neč67, Chap. 3, Thm. 7.6] with [Neč67, Chap. 3, Thm. 7.8]proves m X s =1 k v s k H (Ω) . m +1 X k =1 k N k v k L (Ω) + m X s =1 k v s k L (Ω) . (3.3) Step 3.
This step proves a relation between ( N k ) k =1 ,...,m +1 and sym Curl for a properchoice of v = ( v , . . . , v m ) .Define v = ( v , . . . , v m ) ∈ H (Ω; R m ) by setting for each s ∈ { , . . . , m } the function v s := γ ℓ ,...,ℓ m − with ℓ = · · · = ℓ s − = 1 and ℓ s = · · · = ℓ m − = 2 (with ( ℓ , . . . , ℓ m − ) =(2 , . . . , for s = 1 and ( ℓ , . . . , ℓ m − ) = (1 , . . . , for s = m ). The symmetry of γ proves k Curl γ k L (Ω) . m X s =1 k v s k H (Ω) and m X s =1 k v s k L (Ω) ≈ k γ k L (Ω) . (3.4)With the notation from Step 1 it holds that v s = γ j ( s − = γ j ( s ) and the definition of N k from Step 2 and (3.2) reveal N k +1 v = ( m − k ) /m ( ∂v k +1 /∂x ) − k/m ( ∂v k /∂y ) = (sym Curl γ ) j ( k ) . This leads to m +1 X k =1 k N k v k L (Ω) ≤ k sym Curl γ k L (Ω) . This, (3.4), and an application of (3.3) implies the assertion.Define, for m ≥ , the spaces H (Ω , m −
1) := (cid:8) v ∈ H (Ω; S ( m − (cid:12)(cid:12) ´ Ω v dx = 0 (cid:9) ,Z := { β ∈ H (Ω , m − | sym Curl β = 0 } ,Y := (cid:8) γ ∈ H (Ω , m − (cid:12)(cid:12) ∀ β ∈ Z : (Curl β, Curl γ ) L (Ω) = 0 (cid:9) . (3.5)A computation reveals for m = 2 , that the spaces Z and Y read Z = { γ ∈ H (1) | ∃ c ∈ R , c ∈ R with γ ( x ) = c x + c } ,Y = (cid:8) γ ∈ H (Ω; R ) (cid:12)(cid:12) ´ Ω γ dx = 0 and ´ Ω div γ dx = 0 (cid:9) (3.6)and for m = 3 the space Z reads Z = γ ∈ H (2) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ∃ c , c , c ∈ R , c ∈ R × with γ ( x, y ) = (cid:18) c x + 2 c x c xy + c y + c xc xy + c y + c x c y + 2 c y (cid:19) + c . (3.7)The following theorem generalizes [CGH14, Lemma 3.3] from m = 2 to higher-order ten-sors m > and states that k sym Curl •k L (Ω) defines a norm on Y . Note that k Curl •k L (Ω) = k D •k L (Ω) . 6 heorem 3.2. Any γ ∈ Y satisfies k Curl γ k L (Ω) . k sym Curl γ k L (Ω) . Proof.
Assume for contradiction that the statement does not hold. Then there exists asequence ( γ n ) n ∈ N ∈ Y N with n k sym Curl γ n k L (Ω) ≤ k Curl γ n k L (Ω) = 1 . Since Y ⊆ H ( m − , Poincaré’s inequality implies that all components of γ n are boundedin H (Ω) . Since H (Ω; X ( m − is reflexive and compactly embedded in L (Ω; X ( m − ,there exists a subsequence (not relabelled) with a limit γ ∈ L (Ω; X ( m − , γ n → γ in L (Ω; X ( m − . This and Theorem 3.1 imply k Curl( γ n − γ ℓ ) k L (Ω) . k sym Curl( γ n − γ ℓ ) k L (Ω) + k γ n − γ ℓ k L (Ω) ≤ n + 1 ℓ + k γ n − γ ℓ k L (Ω) → as n, ℓ → ∞ . The Poincaré inequality and the completeness of H (Ω; X ( m − imply the existence of e γ ∈ H (Ω; X ( m − with γ n → e γ in H (Ω; X ( m − and thus γ = e γ . It holds that k sym Curl •k L (Ω) ≤ k Curl •k L (Ω) and, therefore, k sym Curl •k L (Ω) defines a boundedfunctional on H (Ω; X ( m − . Hence, k sym Curl γ k L (Ω) = lim n →∞ k sym Curl γ n k L (Ω) = 0 . (3.8)Let β ∈ Z . Since γ n ∈ Y , the Cauchy inequality reveals (Curl β, Curl γ ) L (Ω) = (Curl β, Curl( γ − γ n )) L (Ω) ≤ k Curl β k L (Ω) k Curl( γ − γ n ) k L (Ω) → as n → ∞ . This and (3.8) lead to γ ∈ Z ∩ Y and therefore γ = 0 . This contradicts k Curl γ k L (Ω) =lim n →∞ k Curl γ n k L (Ω) = 1 and, hence, implies the assertion. Remark . The proof by contradiction from Theorem 3.2does not provide information about the dependency on the domain. A scaling argumentreveals that it does not depend on the size of the domain, but it may depend on its shape.
This section proves a Helmholtz decomposition of L tensors into m th derivatives and thesymmetric part of a Curl in Theorem 4.4. This is a generalization of the Helmholtz decom-position of [BNS07] for fourth-order problems ( m = 2 ). The proof is based on Theorem 4.1below, which characterizes m th-divergence-free smooth functions as symmetric parts ofCurls. Theorem 4.1.
Let m ≥ and τ ∈ C ∞ (Ω; S ( m )) with div m τ = 0 . Then there exists γ ∈ C ∞ (Ω; X ( m − with τ = sym Curl γ. roof. The proof is based on mathematical induction.The base case m = 1 is a classical result [Rud76]. Assume as induction hypothesis thatthe statement holds for ( m − , i.e., for all e τ ∈ C ∞ (Ω; S ( m − with div m − e τ = 0 thereexists γ ∈ C ∞ (Ω; X ( m − with e τ = sym Curl γ .The inductive step is split in five steps. Suppose that τ ∈ C ∞ (Ω; S ( m )) with div m τ = 0 . Step 1.
Then div τ ∈ C ∞ (Ω; X ( m − and div m − div τ = 0 . Let ( j , . . . , j m − ) ∈{ , } m − and σ ∈ S m − . Recall the definition of the divergence from Definition 1. Thesymmetry of τ implies (div τ ) j ,...,j m − = ∂τ j ,...,j m − , /∂x + ∂τ j ,...,j m − , /∂x = ∂τ j σ (1) ,...,j σ ( m − , /∂x + ∂τ j σ (1) ,...,j σ ( m − , /∂x = (div τ ) j σ (1) ,...,j σ ( m − . Hence, div τ ∈ C ∞ (Ω; S ( m − . The induction hypothesis guarantees the existence of β ∈ C ∞ (Ω; X ( m − with div τ = sym Curl β . Step 2.
This step defines some b β ∈ C ∞ (Ω; X ( m )) with div b β = div τ .The definitions of sym and Curl from Section 2 for tensors combine to (sym Curl β ) j ,...,j m − = (card( S m − )) − X σ ∈ S m − ( − j σ ( m − ∂∂x p ( j σ ( m − ) β j σ (1) ,...,j σ ( m − . (4.1)Define b β ∈ C ∞ (Ω; X ( m )) by b β j ,...,j m := ( − p ( j m ) (card( S m − )) − X σ ∈ S m − j σ ( m − = p ( j m ) β j σ (1) ,...,j σ ( m − . (4.2)The definition of b β implies (div b β ) j ,...,j m − = (card( S m − )) − X k =1 ( − p ( k ) ∂∂x k X σ ∈ S m − j σ ( m − = p ( k ) β j σ (1) ,...,j σ ( m − . Since j σ ( m − = p ( k ) if and only if p ( j σ ( m − ) = k , this equals (card( S m − )) − X k =1 X σ ∈ S m − p ( j σ ( m − )= k ( − j σ ( m − ∂∂x p ( j σ ( m − ) β j σ (1) ,...,j σ ( m − = (card( S m − )) − X σ ∈ S m − ( − j σ ( m − ∂∂x p ( j σ ( m − ) β j σ (1) ,...,j σ ( m − and, hence, the combination of the foregoing two displayed formulae with (4.1) leads to div b β = sym Curl β . The combination with Step 1 proves div b β = div τ . Step 3.
Since div( τ − b β ) = 0 , the base case (applied “row-wise” to ( τ − b β ) j ,...,j m − , • )guarantees the existence of γ ∈ C ∞ (Ω; X ( m − with τ − b β = Curl γ . Step 4.
This step shows sym( b β ) = 0 .Let ( j , . . . , j m ) ∈ { , } m be fixed and let N := card( { k ∈ { , . . . , m } | j k = 1 } ) and N := card( { k ∈ { , . . . , m } | j k = 2 } ) be the number of ones and twos. Then M ( j m ) := N − (2 − j m ) and M ( j m ) := N − ( j m − (4.3)8re the numbers of ones and twos in ( j , . . . , j m − ) . Define the index set T := ( ( k , . . . , k m − ) ∈ { , } m − (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) m − X ℓ =1 k ℓ = ( N −
1) + 2( N − ) . This set T contains exactly all indices ( k , . . . , k m − ) with ( N − many ones and ( N − many twos. Note that j σ ( m − = p ( j m ) implies that { j σ ( m − , j m } = { , } and theelements of T are the only indices which appear as indices of β in the sum in (4.2). For j ∈ T , each β j appears M ( j m )! M ( j m )! times in that sum. This and (4.2) yield b β j ,...,j m = ( − p ( j m ) (card( S m − )) − M ( j m )! M ( j m )! X j ∈ T β j . This reveals (sym b β ) j ,...,j m = (card( S m )) − X σ ∈ S m b β j σ (1) ,...,j σ ( m ) = (card( S m )) − (card( S m − )) − × X j ∈ T β j X σ ∈ S m ( − p ( j σ ( m ) ) M ( j σ ( m ) )! M ( j σ ( m ) )! . A reordering of the summands and the definition of M and M in (4.3) leads to X σ ∈ S m ( − p ( j σ ( m ) ) M ( j σ ( m ) )! M ( j σ ( m ) )!= M (1)! M (1)! X σ ∈ S m j σ ( m ) =1 ! − M (2)! M (2)! X σ ∈ S m j σ ( m ) =2 ! = ( N − N !card( { σ ∈ S m | j σ ( m ) = 1 } ) − N ! ( N − { σ ∈ S m | j σ ( m ) = 2 } ) . Since card( { σ ∈ S m | j σ ( m ) = 1 } ) = N card( S m − ) and card( { σ ∈ S m | j σ ( m ) = 2 } ) = N card( S m − ) , this vanishes. This proves sym b β = 0 . Step 5.
Step 4 and τ ∈ C ∞ (Ω; S ( m )) leads to τ = sym( τ ) = sym( τ − b β ) . Step 3 thenyields τ = sym Curl γ and concludes the proof.The following theorem states a Helmholtz decomposition into m th derivatives andsymmetric parts of Curls. The proof uses Theorem 4.1 and a density argument. Thefollowing assumption assumes that the constant in Theorem 3.2 does continuously dependon the domain. To this end, define Ω ε := { x ∈ Ω | dist( x, ∂ Ω) > ε } . (4.4) Assumption 4.2.
There exist sequences ( ε n ) n ∈ N ∈ R N , ( δ n ) n ∈ N ∈ R N , and (Ω ( n ) ) n ∈ N with Ω δ n ⊆ Ω ( n ) ⊆ Ω ε n ⊆ Ω and ε n → and δ n → as n → ∞ , such that the constants C n from Theorem 3.2 with respect to Ω ( n ) are uniformly bounded, sup n ∈ N C n . .Remark . Remark 3.3 implies that Assumption 4.2 is fulfilled on star-shaped domains.Recall the definition of Y from (3.5). 9 heorem 4.4 (Helmholtz decomposition for higher-order derivatives) . If Assumption 4.2is satisfied, then it holds that L (Ω; S ( m )) = D m ( H m (Ω)) ⊕ sym Curl Y and the decomposition is orthogonal in L (Ω; S ( m )) . For any τ ∈ L (Ω; S ( m )) , u ∈ H m (Ω) ,and α ∈ Y with τ = D m u + sym Curl α , the function u ∈ H m (Ω) solves ( D m u, D m v ) = ( τ, D m v ) for all v ∈ H m (Ω) . (4.5) Proof.
Given τ ∈ L (Ω; S ( m )) , let u ∈ H m (Ω) be the solution to (4.5). Define r := τ − D m u ∈ L (Ω; S ( m )) with div m r = 0 .Let ( ε n ) n ∈ N , ( δ n ) n ∈ N , and (Ω ( n ) ) n ∈ N denote the sequences from Assumption 4.2 andlet η n ∈ C ∞ c ( R ) denote the standard mollifier [Eva10] with compact support supp( η n ) in the ball B ε n (0) with radius ε n and centre . Define the regularized function r n := r ∗ η n ∈ C ∞ (Ω; S ( m )) with convolution ∗ . Then r n → r in L (Ω; S ( m )) as n → ∞ . Recallthe definition of Ω ε n from (4.4). Since supp( η n ) ⊆ B ε n (0) and div m r = 0 , it follows (div m r n ) | Ω εn = ( r ∗ D m η n ) | Ω εn = 0 . Since Ω ( n ) ⊆ Ω ε n , Theorem 4.1 guarantees theexistence of γ n ∈ C ∞ (Ω ( n ) ; X ( m − with r n | Ω ( n ) = sym Curl γ n . Recall H ( m − from(3.5) and define Z n := { β n ∈ H (Ω ( n ) , m − | sym Curl β n = 0 } ,Y n := { ζ n ∈ H (Ω ( n ) , m − | ∀ β n ∈ Z n : (Curl β n , Curl ζ n ) L (Ω) = 0 } . Let e γ n ∈ Y n be the orthogonal projection (with respect to (Curl • , Curl • ) L (Ω) ) of γ n to Y n . Then γ n − e γ n ∈ Z n and, hence, sym Curl e γ n = sym Curl γ n = r n | Ω ( n ) . Let ρ n ∈ H (Ω; X ( m − denote the extension of e γ n to Ω with k ρ n k H (Ω) . k e γ n k H (Ω ( n ) ) [LM72,Theorem 8.1]. This, a Poincaré inequality, and Theorem 3.2 together with Assumption 4.2imply k ρ n k H (Ω) . k Curl e γ n k L (Ω ( n ) ) . k sym Curl e γ n k L (Ω ( n ) ) = k r n k L (Ω ( n ) ) . . Since H (Ω; X ( m − is reflexive, there exists a subsequence of ( ρ n ) n ∈ N (again denoted by ρ n ) and γ ∈ H (Ω; X ( m − with ρ n ⇀ γ in H (Ω; X ( m − . Let ϕ ∈ L (Ω; X ( m )) with supp( ϕ ) ⊆ Ω δ n . Since Ω δ n ⊆ Ω ( n ) and therefore sym Curl ρ n | Ω δn = sym Curl e γ n | Ω δn = r n ,it follows ( ϕ, sym Curl γ ) L (Ω) = ( ϕ, r ) L (Ω) + ( ϕ, sym Curl( γ − ρ n )) L (Ω) + ( ϕ, r n − r ) L (Ω) . Since ρ n ⇀ γ in H (Ω; X ( m − and r n → r in L (Ω; S ( m )) and δ n → , this leadsto sym Curl γ = r . Let ρ ∈ Y be the orthogonal projection of γ to Y (with respect to (Curl • , Curl • ) L (Ω) ). Then ρ − γ ∈ Z and, hence, sym Curl ρ = sym Curl γ = r . Thisproves the decomposition.Since Curl is the row-wise application of the standard Curl operator, the L orthog-onality of Curl and ∇ for scalar-valued functions and the symmetry of D m prove the L orthogonality of sym Curl and D m . Subsection 5.1 introduces the weak formulation of problem (1.1) based on the Helmholtzdecomposition from Section 4 and its discretization follows in Subsection 5.2.10 .1 Weak formulation
Recall the definition of the divergence from Section 2 and the definition of Y from (3.5).Let ϕ ∈ H (div m , Ω) with ( − m div m ϕ = f and consider the problem: Seek ( σ, α ) ∈ L (Ω; S ( m )) × Y with ( σ, τ ) L (Ω) + ( τ, sym Curl α ) L (Ω) = ( ϕ, τ ) L (Ω) for all τ ∈ L (Ω; S ( m )) , ( σ, sym Curl β ) L (Ω) = 0 for all β ∈ Y. (5.1)The following theorem states the equivalence of this problem with (1.1). Theorem 5.1 (existence of solutions) . There exists a unique solution ( σ, α ) ∈ L (Ω; S ( m )) × Y to (5.1) with k σ k L (Ω) + k Curl α k L (Ω) . k σ k L (Ω) + k sym Curl α k L (Ω) = k sym ϕ k L (Ω) . (5.2) If Assumption 4.2 is satisfied, then ( σ, α ) satisfies σ = D m u for the solution u ∈ H m (Ω) to (1.1) . Note that σ = D m u is satisfied for any ϕ ∈ H (div m , Ω) with ( − m div m ϕ = f , while sym Curl α = ϕ − D m u depends on the choice of ϕ . Proof of Theorem 5.1.
The inf-sup condition k Curl β k L (Ω) . sup τ ∈ L (Ω; S ( m )) \{ } ( τ, sym Curl β ) L (Ω) k τ k L (Ω) follows from Theorem 3.2. This and Brezzi’s splitting lemma [Bre74] proves the uniqueexistence of a solution to (5.1). Since σ + sym Curl α = sym( ϕ ) , Theorem 3.2 leads to thestability (5.2).If Assumption 4.2 is fulfilled, then the Helmholtz decomposition of Theorem 4.4 holdsand the L orthogonality of σ to sym Curl Y yields the existence of e u ∈ H m (Ω) with σ = D m e u . The orthogonality of Theorem 4.4, ( − m div m ϕ = f , and the symmetry ofthe m th derivative imply for all v ∈ H m (Ω) that ( D m e u, D m v ) L (Ω) = ( ϕ, D m v ) L (Ω) − ( D m v, sym Curl α ) L (Ω) = ( f, v ) . Hence, e u solves (1.1). The discretization of (5.1) employs the discrete spaces X h ( T ) := P k ( T ; S ( m )) ,Y h ( T ) := P k +1 ( T ; X ( m − ∩ Y and seeks σ h ∈ X h ( T ) and α h ∈ Y h ( T ) with ( σ h , τ h ) L (Ω) + ( τ h , sym Curl α h ) L (Ω) = ( ϕ, τ h ) L (Ω) for all τ h ∈ X h ( T ) , ( σ h , sym Curl β h ) L (Ω) = 0 for all β h ∈ Y h ( T ) . (5.3)11 emark . Note that there is no constraint on the polynomial degree k ≥ . A discretiza-tion with the lowest polynomial degree involves only piecewise constant and piecewise affinefunctions for any m ≥ . This should be contrasted to a standard conforming FEM wherethe H m (Ω) conformity causes that the lowest possible polynomial degree is very high (cf.the Argyris FEM with piecewise P functions and 21 local degrees of freedom for m = 2 or the conforming FEM of [Žen70] for arbitrary m with piecewise P m − functions).Discontinuous Galerkin FEMs such as C interior penalty methods [EGH +
02, Bre12] needat least piecewise P functions for m = 2 and piecewise P functions for m = 3 [GN11]. Remark . Since the finite element spaces X h ( T ) and Y h ( T ) differ only in the number ofcomponents and the bilinear forms of (5.3) are similar for all m , an implementation in asingle program which runs for all m is possible. Remark . Since there is no continuity restriction in X h ( T ) be-tween elements, the mass matrix is block diagonal with local mass matrices as sub-blocks.Therefore, the matrix corresponding to the bilinear form ( • , • ) L (Ω) in (5.3) can be directlyinverted. Remark . Problem (5.3) provides an approximation σ h of D m u . If the function u itselfor a lower derivative of u is the quantity of interest, it can be approximated by, e.g., a leastsquares approach. For u h,m := σ h the minimisation of m − X j =0 k u h,j +1 − Du h,j k L (Ω) with respect to ( u h,j ) j =1 ,...,m − over a suitable finite element space results in a series of m Poisson problems and provides an approximation u h, to u . This ansatz can also beemployed to include lower order terms in the system, cf. [Gal15] for a similar approach. Theorem 5.6 (best-approximation result) . There exists a unique solution ( σ h , α h ) ∈ X h ( T ) × Y h ( T ) to (5.3) and it satisfies k σ − σ h k L (Ω) + k sym Curl( α − α h ) k L (Ω) . min τ h ∈ X h ( T ) k σ − τ h k L (Ω) + min β h ∈ Y h ( T ) k sym Curl( α − β h ) k L (Ω) . (5.4)If the solution is sufficiently smooth, say σ ∈ H k +1 (Ω; S (2)) and α ∈ H k +2 (Ω; R ) , thisyields a convergence rate of O ( h k +1 ) . Remark ϕ ) . Given a right-hand side f , the discretization (5.3) requiresthe knowledge of a function ϕ ∈ H (div m , Ω) with ( − m div m ϕ = f . This can be computedby an integration of f – manually for a simple f or numerically for a more complicated f . This can be done in parallel. However, the numerical experiments of Section 7 and thebest-approximation result in Theorem 5.6 suggest that the magnitude of the error heavilydepends on the choice of ϕ (which determines sym Curl α ). In Section 7, the error can bedrastically reduced by defining ϕ by ϕ = ∆ − ∇ ∆ − ∇ ∆ − f and approximate ∆ − withstandard finite elements (see Section 7 for more details). Proof of Theorem 5.6.
Since sym Curl Y h ( T ) ⊆ X h ( T ) , Theorem 3.2 proves the inf-supcondition k Curl β h k L (Ω) . sup τ h ∈ X h ( T ) \{ } ( τ h , sym Curl β h ) L (Ω) k τ h k L (Ω) for all β h ∈ Y h ( T ) . Brezzi’s splitting lemma [Bre74] therefore leads to the unique existence of a solution ofproblem (5.3). This, the conformity of the discretization, and standard arguments formixed FEMs [BBF13] lead to the best-approximation result (5.4).12efine the space of discrete orthogonal derivatives as W h ( T ) := { τ h ∈ X h ( T ) | ∀ β h ∈ Y h ( T ) : ( τ h , sym Curl β h ) L (Ω) = 0 } . (5.5)The following lemma proves a projection property. Lemma 5.8 (projection property) . Let τ ∈ L (Ω; S ( m )) with ( τ, sym Curl β ) L (Ω) = 0 for all β ∈ Y. Then Π X h ( T ) τ ∈ W h ( T ) . If T ⋆ is an admissible refinement of T and τ ⋆ ∈ W h ( T ⋆ ) , then Π X h ( T ) τ ⋆ ∈ W h ( T ) .Proof. Let β h ∈ Y h ( T ) . Since sym Curl β h ∈ X h ( T ) , the conformity Y h ( T ) ⊆ Y implies (Π X h ( T ) τ, sym Curl β h ) = ( τ, sym Curl β h ) = 0 . The same arguments apply to τ ⋆ ∈ W h ( T ⋆ ) . For m = 2 , problem (1.1) becomes the biharmonic problem ∆ u = f . This problem arisesin the theory of Kirchhoff plates with clamped boundary. In this situation, the Helmholtzdecomposition of Theorem 4.4 is already proved in [BNS07].The discrete spaces in (5.3) for m = 2 read X h = P k ( T ; S (2)) with S (2) the space ofsymmetric × matrices and Y h = P k +1 ( T ; R ) ∩ Y with Y defined in (3.6). For platebending problems, [Mor68] introduced a P non-conforming finite element method (alsocalled Morley FEM) with non-conforming finite element space V M ( T ) := v h ∈ P ( T ) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) v h is continuous at the interior nodes and vanishes atboundary nodes; ∇ NC v h is continuous at the interioredges’ midpoints and vanishes at the midpoints ofboundary edges . The discrete Helmholtz decomposition [CGH14] P ( T ; S (2)) = D V M ( T ) ⊕ sym Curl (cid:0) P ( T ; R ) ∩ Y (cid:1) . shows for k = 0 the relation D V M ( T ) = W h ( T ) with W h ( T ) from (5.5) and, hence, thesolution σ h to (5.3) is a piecewise Hessian of a Morley function. If ϕ satisfies div ϕ = f also in the dual space of V M ( T ) , then the solution σ h ∈ X h ( T ) of (5.3) coincides with thepiecewise Hessian of the solution of the Morley FEM.For m = 3 , problem (1.1) becomes the triharmonic problem − ∆ u = f . Sixth-orderequations arise in the description of the motion of thin viscous droplets [BLN04] or of theoxidation of silicon in superconductor devices [Kin89]. For the triharmonic problem, thediscrete spaces read X h = P k ( T ; S (3)) and Y h = P k +1 ( T ; R × ) ∩ Y with Y defined in (3.5).The orthogonality onto Z implied by the definition of Y can be implemented by Lagrangemultipliers and with the knowledge of Z from (3.7). This section defines the adaptive algorithm and proves its quasi-optimal convergence.13 .1 Adaptive algorithm and optimal convergence rates
Let T denote some initial shape-regular triangulation of Ω , such that each triangle T ∈ T is equipped with a refinement edge E T ∈ E ( T ) . We assume that T fulfils the followinginitial condition. Definition 2 (initial condition) . All
T, K ∈ T with T ∩ K = E ∈ E and with refinementedges E T ∈ E ( T ) and E K ∈ E ( K ) satisfy: If E T = E , then E K = E T . If E K = E , then E T = E K .Given an initial triangulation T , the set of admissible triangulations T is defined as theset of all regular triangulations which can be created from T by newest-vertex bisection(NVB) [Ste08]. Let T ( N ) denote the subset of all admissible triangulations with at most card( T ) + N triangles. The adaptive algorithm involves the overlay of two admissibletriangulations T , T ⋆ ∈ T , which reads T ⊗ T ⋆ := { T ∈ T ∪ T ⋆ | ∃ K ∈ T , K ⋆ ∈ T ⋆ with T ⊆ K ∩ K ⋆ } . (6.1)The adaptive algorithm is based on separate marking. Given a triangulation T ℓ , definefor all T ∈ T ℓ the local error estimator contributions by λ ( T ℓ , T ) := k h T curl NC σ h k L ( T ) + h T X E ∈ E ( T ) k [ σ h ] E · τ E k L ( E ) ,µ ( T ) := k sym( ϕ ) − Π k sym( ϕ ) k L ( T ) and the global error estimators by λ ℓ := λ ( T ℓ , T ℓ ) with λ ( T ℓ , M ) := X T ∈ M λ ( T ℓ , T ) for all M ⊆ T ℓ ,µ ℓ := µ ( T ℓ ) with µ ( M ) := X T ∈ M µ ( T ) for all M ⊆ T ℓ . The adaptive algorithm is driven by these two error estimators and runs the following loop.
Algorithm 6.1 (AFEM for higher-order problems) . Input:
Initial triangulation T , pa-rameters < θ A ≤ , < ρ B < , < κ . for ℓ = 0 , , , . . . do Solve.
Compute solution ( σ ℓ , α ℓ ) ∈ X h ( T ℓ ) × Y h ( T ℓ ) of (5.3) with respect to T ℓ . Estimate.
Compute estimator contributions (cid:0) λ ( T ℓ , T ) (cid:1) T ∈ T ℓ and (cid:0) µ ( T ) (cid:1) T ∈ T ℓ . if µ ℓ ≤ κλ ℓ then Mark.
The Dörfler marking chooses a minimal subset M ℓ ⊆ T ℓ such that θ A λ ℓ ≤ λ ℓ ( T ℓ , M ℓ ) Refine.
Generate the smallest admissible refinement T ℓ +1 of T ℓ in whichat least all triangles in M ℓ are refined. else Mark.
Compute an admissible triangulation T ∈ T with µ T ≤ ρ B µ ℓ . Refine.
Generate the overlay T ℓ +1 of T ℓ and T (cf. (6.1)). end ifend forOutput: Sequence of triangulations ( T ℓ ) ℓ ∈ N and discrete solutions ( σ ℓ , α ℓ ) ℓ ∈ N . (cid:7) µ ℓ > κλ ℓ can be realized by the algorithm Approx from [CR15, BDD04], i.e. the threshold second algorithm [BD04] followed by a completionalgorithm.For s > and ( σ, α, ϕ ) ∈ L (Ω; S ( m )) × Y × H (div m , Ω) , define | ( σ, α, ϕ ) | A s := sup N ∈ N N s inf T ∈ T ( N ) (cid:16) (cid:13)(cid:13) σ − Π X h ( T ) σ (cid:13)(cid:13) L (Ω) + inf β T ∈ Y h ( T ) k sym Curl( α − β T ) k L (Ω) + k ϕ − Π X h ( T ) ϕ k L (Ω) (cid:17) . Remark . A “row-wise” application of [Vee14, Theo-rem 3.2] proves | ( σ, α, ϕ ) | A s ≈ | ( σ, α, ϕ ) | A ′ s := sup N ∈ N N s inf T ∈ T ( N ) (cid:16) k σ − Π X h ( T ) σ k L (Ω) + (cid:13)(cid:13) sym(Curl α ) − Π X h ( T ) sym(Curl α ) (cid:13)(cid:13) L (Ω) + (cid:13)(cid:13) sym( ϕ ) − Π X h ( T ) sym( ϕ ) (cid:13)(cid:13) L (Ω) (cid:17) . In the following, we assume that the following axiom (B1) holds for the algorithm usedin the step
Mark for µ ℓ > κλ ℓ . For the algorithm Approx , this assumption is a consequenceof Axioms (B2) and (SA) from Subsection 6.5 [CR15].
Assumption 6.3 ((B1) optimal data approximation) . Assume that | ( σ, α, ϕ ) | A s is finite.Given a tolerance Tol , the algorithm used in Mark in the second case ( µ ℓ > κλ ℓ ) inAlgorithm 6.1 computes T ⋆ ∈ T with card( T ⋆ ) − card( T ) . Tol − / (2 s ) and µ ( T ⋆ ) ≤ Tol . The following theorem states optimal convergence rates of Algorithm 6.1.
Theorem 6.4 (optimal convergence rates of AFEM) . For < ρ B < and sufficientlysmall < κ and < θ < , Algorithm 6.1 computes sequences of triangulations ( T ℓ ) ℓ ∈ N and discrete solutions ( σ ℓ , α ℓ ) ℓ ∈ N for the right-hand side ϕ of optimal rate of convergencein the sense that (card( T ℓ ) − card( T )) s (cid:16) k σ − σ ℓ k L (Ω) + k sym Curl( α − α ℓ ) k L (Ω) (cid:17) . | ( σ, α, ϕ ) | A s . The proof follows from the abstract framework of [CR15], under the assumptions(A1)–(A4), which are proved in Subsections 6.2–6.4, the assumption (B1), which followsfrom (B2) and (SA) from Subection 6.5 below for the algorithm
Approx , and efficiency of p λ + µ , which follows from the standard bubble function technique of [Ver96]. The following two theorems follow from the structure of the error estimator λ . Theorem 6.5 (stability) . Let T ⋆ be an admissible refinement of T and M ⊆ T ∩ T ⋆ . Let ( σ T ⋆ , α T ⋆ ) ∈ X h ( T ⋆ ) × Y h ( T ⋆ ) and ( σ T , α T ) ∈ X h ( T ) × Y h ( T ) be the respective discretesolutions to (5.3) . Then, | λ ( T ⋆ , M ) − λ ( T , M ) | . k σ T ⋆ − σ T k L (Ω) . roof. This follows with triangle inequalities, inverse inequalities and the trace inequalityfrom [BS08, p. 282] as in [CKNS08, Proposition 3.3].
Theorem 6.6 (reduction) . Let T ⋆ be an admissible refinement of T . Then there exists < ρ < and Λ < ∞ such that λ ( T ⋆ , T ⋆ \ T ) ≤ ρ λ ( T , T \ T ⋆ ) + Λ k σ T ⋆ − σ T k . Proof.
This follows with a triangle inequality and the mesh-size reduction property h T ⋆ | T ≤ h T | T / for all T ∈ T ⋆ \ T as in [CKNS08, Corollary 3.4]. The following theorem proves discrete reliability, i.e., the difference between two discretesolutions is bounded by the error estimators on refined triangles only.
Theorem 6.7 (discrete reliability) . Let T ⋆ be an admissible refinement of T with respectivediscrete solutions ( σ T ⋆ , α T ⋆ ) ∈ X h ( T ⋆ ) × Y h ( T ⋆ ) and ( σ T , α T ) ∈ X h ( T ) × Y h ( T ) of (5.3) .Then, k σ T − σ T ⋆ k + k sym Curl( α T − α T ⋆ ) k L (Ω) . λ ( T , T \ T ⋆ ) + µ ( T , T \ T ⋆ ) . Proof.
Recall the definition of W h ( T ⋆ ) from (5.5). Since σ T − σ T ⋆ ∈ X h ( T ⋆ ) , there exist p T ⋆ ∈ W h ( T ⋆ ) and r T ⋆ ∈ Y h ( T ⋆ ) with σ T − σ T ⋆ = p T ⋆ + sym Curl r T ⋆ . The discrete error canbe split as k σ T − σ T ⋆ k L (Ω) = ( σ T − σ T ⋆ , p T ⋆ ) L (Ω) + ( σ T − σ T ⋆ , sym Curl r T ⋆ ) L (Ω) . (6.2)The projection property, Lemma 5.8, proves Π X h ( T ) p T ⋆ ∈ W h ( T ) . Hence, problem (5.3)implies that the first term of the right-hand side equals ( σ T − σ T ⋆ , p T ⋆ ) L (Ω) = (Π X h ( T ) ϕ − ϕ, p T ⋆ ) L (Ω) = (Π X h ( T ) ϕ − Π X h ( T ⋆ ) ϕ, p T ⋆ ) L (Ω) . For any triangle T ∈ T ∩ T ⋆ , it holds (Π X h ( T ) ϕ − Π X h ( T ⋆ ) ϕ ) | T = 0 . Since T ⋆ is a refinementof T , this implies (Π X h ( T ) ϕ − Π X h ( T ⋆ ) ϕ, p T ⋆ ) L (Ω) ≤ k Π X h ( T ) ϕ − Π X h ( T ⋆ ) ϕ k T \ T ⋆ k p T ⋆ k L (Ω) ≤ k ϕ − Π X h ( T ) ϕ k T \ T ⋆ k p T ⋆ k L (Ω) . Let r T ∈ Y h ( T ) denote the quasi interpolant from [SZ90] of r T ⋆ which satisfies theapproximation and stability properties k h − T ( r T ⋆ − r T ) k L (Ω) + k D ( r T ⋆ − r T ) k L (Ω) . k Dr T ⋆ k L (Ω) and ( r T ) | E = ( r T ⋆ ) | E for all edges E ∈ E ( T ) ∩ E ( T ⋆ ) . Since σ T ∈ W h ( T ) and σ T ⋆ ∈ W h ( T ⋆ ) ,the symmetry of σ T implies ( σ T − σ T ⋆ , sym Curl r T ⋆ ) L (Ω) = ( σ T , sym Curl( r T ⋆ − r T )) L (Ω) = ( σ T , Curl( r T ⋆ − r T )) L (Ω) . (6.3)An integration by parts leads to ( σ T , Curl( r T ⋆ − r T )) L (Ω) = − (curl NC σ T , r T ⋆ − r T ) L (Ω) + X E ∈ E ( T ) ˆ E [ σ T · τ E ] E ( r T ⋆ − r T ) ds. T ∈ T ∩ T ⋆ , any edge E ∈ E ( T ) satisfies E ∈ E ( T ) ∩ E ( T ⋆ ) and, hence, ( r T ) | T = ( r T ⋆ ) | T for all T ∈ T ∩ T ⋆ . This, the Cauchy inequality, the approximation andstability properties of the quasi interpolant, and the trace inequality from [BS08, p. 282]lead to − (curl NC σ T , r T ⋆ − r T ) L (Ω) + X E ∈ E ˆ E [ σ T · τ E ] E ( r T ⋆ − r T ) ds . (cid:18) k h T curl NC σ T k T \ T ⋆ + s X E ∈ E ( T ) \ E ( T ⋆ ) h T k [ σ T · τ E ] E k L ( E ) (cid:19) k Curl r T ⋆ k L (Ω) . (6.4)The combination of the previous displayed inequalities yields k σ T − σ T ⋆ k L (Ω) . λ ( T , T \ T ⋆ ) + µ ( T , T \ T ⋆ ) . Since
Curl α T = Π X h ( T ) ϕ − σ T and Curl α T ⋆ = Π X h ( T ⋆ ) ϕ − σ T ⋆ , the triangle inequality yieldsthe assertion. Remark . The convergence of σ T ⋆ and α T ⋆ , whichis a consequence of the a priori error estimate of Theorem 5.6, and the discrete reliabilityof Theorem 6.7 imply the reliability k σ − σ T k L (Ω) + k sym Curl( α − α T ) k L (Ω) . λ ℓ + µ ℓ . The following theorem proves quasi-orthogonality of the discretization (5.3).
Theorem 6.9 (general quasi-orthogonality) . Let ( T j | j ∈ N ) be some sequence of tri-angulations with discrete solutions ( σ j , α j ) ∈ X h ( T j ) × Y h ( T j ) to (5.3) and let ℓ ∈ N .Then, ∞ X j = ℓ (cid:16) k σ j − σ j − k + k sym Curl( α j − α j − ) k (cid:17) . λ ℓ − + µ ℓ − . Proof.
The projection property, Lemma 5.8, proves Π X h ( T j − ) σ j ∈ W h ( T j − ) with W h ( T j − ) from (5.5). Hence, problem (5.3) leads to ( σ j − , σ j − σ j − ) L (Ω) = ( ϕ, Π X h ( T j − ) σ j − σ j − ) L (Ω) , ( σ j , σ j − σ j − ) L (Ω) = ( ϕ, σ j ) − ( ϕ, Π X h ( T j − ) σ j ) L (Ω) . The subtraction of these two equations and an index shift leads, for any M ∈ N with M > ℓ , to M X j = ℓ k σ j − σ j − k L (Ω) = M X j = ℓ ( ϕ, σ j − Π X h ( T j − ) σ j ) L (Ω) − M X j = ℓ ( ϕ, Π X h ( T j − ) σ j ) L (Ω) + M − X j = ℓ − ( ϕ, σ j ) L (Ω) = ( ϕ, σ ℓ − − σ M ) L (Ω) + 2 M X j = ℓ ( ϕ, σ j − Π X h ( T j − ) σ j ) L (Ω) . (6.5)17ince σ j − Π X h ( T j − ) σ j ∈ X h ( T j ) is L -orthogonal to X h ( T j − ) , a Cauchy and a weightedYoung inequality imply M X j = ℓ ( ϕ, σ j − Π X h ( T j − ) σ j ) L (Ω) = 2 M X j = ℓ (Π X h ( T j ) ϕ − Π X h ( T j − ) ϕ, σ j − Π X h ( T j − ) σ j ) L (Ω) ≤ M X j = ℓ k Π X h ( T j ) ϕ − Π X h ( T j − ) ϕ k L (Ω) + 12 M X j = ℓ k σ j − Π X h ( T j − ) σ j k L (Ω) . (6.6)The orthogonality Π X h ( T j ) ϕ − Π X h ( T j − m ) ϕ ⊥ L (Ω) X h ( T j − m ) for all ≤ m ≤ j and thedefinition of µ ℓ proves M X j = ℓ k Π X h ( T j ) ϕ − Π X h ( T j − ) ϕ k L (Ω) = k Π X h ( T M ) ϕ − Π X h ( T ℓ − ) ϕ k L (Ω) = k Π X h ( T M ) ( ϕ − Π X h ( T ℓ − ) ϕ ) k L (Ω) ≤ µ ℓ − . (6.7)The combination of (6.5)–(6.7) and k σ j − Π X h ( T j − ) σ j k L (Ω) ≤ k σ j − σ j − k L (Ω) leads to M X j = ℓ k σ j − σ j − k L (Ω) ≤ µ ℓ − + ( ϕ, σ ℓ − − σ M ) L (Ω) . (6.8)The arguments of (6.3)–(6.4) prove (sym Curl( α M − α ℓ − ) , σ ℓ − ) L (Ω) . λ ℓ − k Curl( α M − α ℓ − ) k L (Ω) . The discrete problem (5.3), the discrete reliability k sym Curl( α M − α ℓ − ) k L (Ω) . λ ℓ − + µ ℓ − from Theorem 6.7, and Theorem 3.2 therefore lead to ( σ ℓ − − σ M , Π X h ( T ℓ − ) ϕ ) L (Ω) = ( σ ℓ − − σ M , σ ℓ − + Curl α ℓ − ) L (Ω) = ( σ ℓ − − σ M , σ ℓ − ) L (Ω) = (sym Curl( α M − α ℓ − ) , σ ℓ − ) L (Ω) . λ ℓ − k Curl( α M − α ℓ − ) k L (Ω) . λ ℓ − + µ ℓ − . (6.9)This and a further application of Theorem 6.7 leads to ( ϕ, σ ℓ − − σ M ) L (Ω) = ( ϕ − Π X h ( T ℓ − ) ϕ, σ ℓ − − σ M ) L (Ω) + ( σ ℓ − − σ M , Π X h ( T ℓ − ) ϕ ) L (Ω) . k ϕ − Π X h ( T ℓ − ) ϕ k L (Ω) k σ ℓ − − σ M k L (Ω) + ( λ ℓ − + µ ℓ − ) L (Ω) . λ ℓ − + µ ℓ − . (6.10)The combination of (6.8) with (6.10) implies M X j = ℓ k σ j − σ j − k L (Ω) . λ ℓ − + µ ℓ − . (6.11)18he Young inequality, the triangle inequality, and sym Curl α j = Π X h ( T j ) ϕ − σ j imply M X j = ℓ k sym Curl( α j − α j − ) k L (Ω) ≤ M X j = ℓ k σ j − σ j − k L (Ω) + 2 M X j = ℓ k Π X h ( T j ) ϕ − Π X h ( T j − ) ϕ k L (Ω) . Since
M > ℓ is arbitrary, the combination with (6.7) and (6.11) yields the assertion.
The following theorem states quasi-monotonicity and sub-additivity for the data-approx-imation error estimator µ . This theorem implies that Assumption 6.3 is satisfied if thealgorithm Approx from [BD04, BDD04, CR15] is used in the second marking step ( µ ℓ ≥ κλ ℓ ) in Algorithm 6.1 [CR15]. Theorem 6.10 ((B2) quasi-monotonicity and (SA) sub-additivity) . Any admissible re-finement T ⋆ of T satisfies µ ( T ⋆ ) ≤ µ ( T ) and X T ∈ T ⋆ T ⊆ K µ ( T ) ≤ µ ( K ) for all K ∈ T . Proof.
This follows directly from the definition of µ . This section is devoted to numerical experiments for the plate problem ∆ u = f and thesixth-order problem − ∆ u = f . The discretization (5.3) is realized for k = 0 , for theplate problem and for k = 0 , , for the sixth-order problem. The experiments comparethe errors and error estimators on a sequence of uniformly red-refined triangulations (thatis, the midpoints of the edges of a triangle are connected; this generates four new tri-angles) with the errors and error estimators on a sequence of triangulations created byAlgorithm 6.1 with bulk parameter θ = 0 . and κ = 0 . and ρ = 0 . .The convergence history plots are logarithmically scaled and display the error k σ − σ h k L (Ω) against the number of degrees of freedom (ndof) of the linear system resultingfrom the Schur complement. m = 2 The exact solution to ∆ u ( x, y ) = f ( x, y ) :=24( x − x + x + y − y + y )+ 2(2 − x + 12 x )(2 − y + 12 y ) with clamped boundary conditions u | ∂ Ω = ( ∂u/∂ν ) | ∂ Ω = 0 reads u ( x, y ) = x (1 − x ) y (1 − y ) . − − − − − − −
11 1 0 . ndof η for k = 0 , uniform k σ − σ h k L for k = 0 , uniform η for k = 0 , adaptive k σ − σ h k L for k = 0 , adaptive η for k = 1 , uniform k σ − σ h k L for k = 1 , uniform η for k = 1 , adaptive k σ − σ h k L for k = 1 , adaptive Figure 2: Errors and error estimators for the experiment on the square from Subsection 7.1.
Define ϕ = ( ϕ jk ) ≤ j,k ≤ ∈ H (div , Ω) by ϕ := 24( x / − x /
10 + x /
30) + ( x − x + x )(2 − y + 12 y ) ,ϕ := 24( y / − y /
10 + y /
30) + ( y − y + y )(2 − x + 12 x ) ,ϕ := ϕ := 0 . Then div ϕ = f and ϕ is an admissible right-hand side for (5.3).The errors k σ − σ h k L (Ω) and error estimators p λ + µ are plotted in Figure 2 versusthe degrees of freedom. The errors and error estimators show an equivalent behaviourwith an overestimation factor of approximately 10. The errors and error estimators showa convergence rate of ndof − / for k = 0 and of ndof − for k = 1 on the sequence ofuniformly red-refined triangulations as well as on the sequence of triangulations generatedby Algorithm 6.1. All marking steps in Algorithm 6.1 for k = 0 , applied the Dörflermarking ( µ ℓ ≤ κλ ℓ ). m = 2 This subsection considers the problem ∆ u = 1 on the L-shaped domain Ω := ( − , \ ([0 , × [ − , with clamped boundary conditions u | ∂ Ω = ( ∂u/∂ν ) | ∂ Ω = 0 and unknown solution. Define the right-hand side ϕ ∈ H (div , Ω) with div ϕ = 1 by ϕ ( x, y ) := (cid:18) x / y / (cid:19) . The error estimators p λ + µ are plotted in Figure 3 versus the degrees of freedom.For uniform mesh-refinement the convergence rate of the error estimator for k = 1 is ndof − / . The convergence rate for k = 0 is slightly larger, but the size of the errorestimator is larger than for k = 1 . This suggests that the observed higher convergence rateis a preasymptotic effect. On the sequences of triangulations generated by Algorithm 6.1,20 − − − − − −
11 1 0 . ndof η for k = 0 , uniform η for k = 0 , adaptive η for k = 1 , uniform η for k = 1 , adaptive Figure 3: Error estimators for the experiment on the L-shaped domain from Subsection 7.2.Figure 4: Adaptively refined triangulations for k = 0 with 1096 nodes (2195 dofs) and for k = 1 with 1077nodes (5114 dofs) for the experiment on the L-shaped domain from Subsection 7.2. − − − − − − − . . . ndof η for k = 0 , uniform η for k = 0 , adaptive η for k = 1 , uniform η for k = 1 , adaptive η for k = 2 , uniform η for k = 2 , adaptive Figure 5: Errors and error estimators for the experiment on the square for m = 3 from Subsection 7.3.The dashed lines correspond to the right-hand side generated by the solution of three successive Poissonproblems. the error estimators show the optimal convergence rates of ndof − / and ndof − for k = 0 and k = 1 , respectively. Figure 4 displays triangulations with approximately 1000 verticesgenerated by Algorithm 6.1 for k = 0 and k = 1 . A stronger refinement towards there-entrant corner is clearly visible. The marking with respect to the data-approximation( µ ℓ > κλ ℓ in Algorithm 6.1) is only applied at the first two levels for k = 0 . All othermarking steps for k = 0 , use the Dörfler marking ( µ ℓ ≤ κλ ℓ ). m = 3 In this subsection, let
Ω = (0 , be the unit square and u ∈ H (Ω) be defined by u ( x, y ) = x (1 − x ) y (1 − y ) with corresponding right-hand side f := − ∆ u . Let ϕ = ( ϕ jkℓ ) ≤ j,k,ℓ ≤ ∈ H (div , Ω) bedefined by ϕ ( x, y ) := − ˆ x ˆ s ˆ t f ( ξ, y ) dξ dt ds,ϕ ( x, y ) := − ˆ y ˆ s ˆ t f ( x, ξ ) dξ dt ds,ϕ := ϕ := ϕ := ϕ := ϕ := ϕ := 0 . (7.1)Then − div ϕ = f and ϕ is an admissible right-hand side for (5.3).The errors k σ − σ h k L (Ω) and error estimators p λ + µ are plotted in Figure 5 versusthe number of degrees of freedom. The errors show the optimal convergence rates of ndof − / , ndof − , and ndof − / for k = 0 , , for uniform refinement as well as for thesequence of triangulations generated by Algorithm 6.1. The error estimators for k = 0 , , show an equivalent behaviour as the respective errors with an overestimation between 3and 9.Although the convergence rates are optimal, one has to consider that the H -seminormof the exact solution k σ k L (Ω) is approximately × − . That means that the relativeerrors for k = 1 (resp. k = 2 ) are larger than up to (resp. ) degrees of22 igure 6: Adaptively refined triangulations for k = 0 with 1460 nodes (4386 dofs), for k = 1 with 1555nodes (19369 dofs), and for k = 2 with 1547 nodes (40833 dofs) for the experiment on the square fromSubsection 7.3. freedom and for k = 0 , they do not even reach this threshold. While the L norm of thefunction σ of interest is approximately − , the L norm of ϕ (and thus k Curl α k L (Ω) ) isapproximately 80. The best-approximation result (5.4) therefore seems to suffer from thelarge term inf β h ∈ Y h ( T ) k Curl( α − β h ) k L (Ω) on the right-hand side.A second choice for the right-hand side ϕ should indicate one possibility to decrease theerror. To this end, define e ϕ := ∇ w with ( w , w , w ) ∈ H (Ω) × H (Ω; R ) × H (Ω; R × ) the solution of ( ∇ w , ∇ v ) L (Ω) = ( f, v ) L (Ω) for all v ∈ H (Ω)( ∇ w , ∇ v ) L (Ω) = ( ∇ w , v ) L (Ω) for all v ∈ H (Ω; R )( ∇ w , ∇ v ) L (Ω) = ( ∇ w , v ) L (Ω) for all v ∈ H (Ω; R × ) . (7.2)Then − div e ϕ = f and the computations are performed with the approximation e ϕ h of e ϕ computed by the approximation of the Poisson problems (7.2) by standard conformingFEMs of degree k . The errors for this right-hand side are included in Figure 5 for k = 0 , , with dashed lines. The errors show the optimal convergence rates and the size of the errorsare reduced by a factor between and compared to the errors for the right-hand sidegiven by (7.1). In this situation, the error is below for all triangulations.Figure 6 displays triangulations with approximately 1500 vertices generated by Algo-rithm 6.1 for k = 0 , , . Although the solution is smooth, a strong refinement towards thecorner (1 , can be observed for k = 0 . For k = 1 , there is a slight refinement towards thecorner (1 , , while for k = 2 , the refinement is nearly uniform. Since the relative errorsfor k = 0 , are still over on these triangulations, the discrete solution probably donot reflect the behaviour of the exact smooth solution. However, the convergence rates areoptimal and the error is slightly smaller compared with the uniform refinement. This is inagreement with Theorem 6.4.All marking steps in Algorithm 6.1 for k = 0 , , used the Dörfler marking ( µ ℓ ≤ κλ ℓ ). m = 3 This section considers the problem: Find u ∈ H (Ω) with − ∆ u = 1 − − − − − − − . . . ndof η for k = 0 , uniform η for k = 0 , adaptive η for k = 1 , uniform η for k = 1 , adaptive η for k = 2 , uniform η for k = 2 , adaptive Figure 7: Errors and error estimators for the experiment on the L-shaped domain for m = 3 from Subsec-tion 7.4.Figure 8: Adaptively refined triangulations for k = 0 with 1118 nodes (3360 dofs), for k = 1 with 1005nodes (11679 dofs), and for k = 2 with 1004 nodes (25875 dofs) for the experiment on the L-shaped domainfrom Subsection 7.4. and homogeneous Dirichlet boundary conditions on the L-shaped domain Ω := ( − , \ ([0 , × [ − , . Let ϕ = ( ϕ jkℓ ) ≤ j,k,ℓ ≤ ∈ H (div , Ω) be defined by ϕ ( x, y ) := − x / ,ϕ := − y / ,ϕ := ϕ := ϕ := ϕ := ϕ := ϕ := 0 Then − div ϕ = 1 and ϕ is an admissible right-hand side for (5.3).Since the exact solution is not known, only the error estimators p λ + µ are plottedin Figure 7 for k = 0 , , on a sequence of uniformly red-refined triangulations and ona sequence generated by Algorithm 6.1. On the sequence of uniformly refined meshes,the error estimators for k = 1 , show a convergence rate of ndof − / , while the errorestimator for k = 0 converges with rate / . However, this error estimator is of larger sizethan the error estimators for k = 1 , and it is therefore expected that the higher rate isa preasymptotic effect. Algorithm 6.1 leads to the optimal convergence rates of ndof − / for k = 0 , ndof − for k = 1 , and ndof − / for k = 2 .Figure 8 displays triangulations with approximately 1000 vertices generated by Algo-rithm 6.1 for k = 0 , , . The strong refinement towards the re-entrant corner is clearly24isible for k = 1 , , while for k = 0 the refinement is quasi-uniform. This is in agreementwith the observed convergence rate for k = 0 and the interpretation that the behaviour ofthe exact solution is not reflected in the discrete solution up to this number of degrees offreedom. The marking with respect to the data-approximation ( µ ℓ > κλ ℓ in Algorithm 6.1)is only applied at levels 1 and 2 for k = 0 . All other marking steps for k = 0 , , use theDörfler marking ( µ ℓ ≤ κλ ℓ ). References [AF89] D. N. Arnold and R. S. Falk. A uniformly accurate finite element method forthe Reissner-Mindlin plate.
SIAM J. Numer. Anal. , 26(6):1276–1290, 1989.[BBF13] D. Boffi, F. Brezzi, and M. Fortin.
Mixed Finite Element Methods and Applica-tions , volume 44 of
Springer Series in Computational Mathematics . Springer,Heidelberg, 2013.[BD04] P. Binev and R. DeVore. Fast computation in adaptive tree approximation.
Numer. Math. , 97(2):193–217, 2004.[BDD04] P. Binev, W. Dahmen, and R. DeVore. Adaptive finite element methods withconvergence rates.
Numer. Math. , 97(2):219–268, 2004.[BLN04] J. W. Barrett, S. Langdon, and R. Nürnberg. Finite element approximation of asixth order nonlinear degenerate parabolic equation.
Numer. Math. , 96(3):401–434, 2004.[BNS07] L. Beirão da Veiga, J. Niiranen, and R. Stenberg. A posteriori error estimatesfor the Morley plate bending element.
Numer. Math. , 106(2):165–179, 2007.[Bre74] F. Brezzi. On the existence, uniqueness and approximation of saddle-pointproblems arising from Lagrangian multipliers.
Rev. Française Automat. Infor-mat. Recherche Opérationnelle Sér. Rouge , 8(R-2):129–151, 1974.[Bre12] S. C. Brenner. C interior penalty methods. In Frontiers in NumericalAnalysis—Durham 2010 , volume 85 of
Lect. Notes Comput. Sci. Eng. , pages79–147. Springer, Heidelberg, 2012.[BS08] S. C. Brenner and L. R. Scott.
The Mathematical Theory of Finite ElementMethods , volume 15 of
Texts in Applied Mathematics . Springer Verlag, NewYork, Berlin, Heidelberg, 3 edition, 2008.[CGH14] C. Carstensen, D. Gallistl, and J. Hu. A discrete Helmholtz decomposition withMorley finite element functions and the optimality of adaptive finite elementschemes.
Comput. Math. Appl. , 68(12):2167–2181, 2014.[Cia78] Ph. G. Ciarlet.
The Finite Element Method for Elliptic Problems . Studiesin Mathematics and its Applications, Vol. 4. North-Holland Publishing Co.,Amsterdam-New York-Oxford, 1978.[CKNS08] J. M. Cascon, Ch. Kreuzer, R. H. Nochetto, and K. G. Siebert. Quasi-optimalconvergence rate for an adaptive finite element method.
SIAM J. Numer. Anal. ,46(5):2524–2550, 2008. 25CR73] M. Crouzeix and P.-A. Raviart. Conforming and nonconforming finite elementmethods for solving the stationary Stokes equations. I.
Rev. Française Automat.Informat. Recherche Opérationnelle Sér. Rouge , 7(R-3):33–75, 1973.[CR15] C. Carstensen and H. Rabus. Axioms of adaptivity for separate marking.preprint, arXiv:1606.02165, 2016.[EGH +
02] G. Engel, K. Garikipati, T. J. R. Hughes, M. G. Larson, L. Mazzei, and R. L.Taylor. Continuous/discontinuous finite element approximations of fourth-orderelliptic problems in structural and continuum mechanics with applications tothin beams and plates, and strain gradient elasticity.
Comput. Methods Appl.Mech. Engrg. , 191(34):3669–3750, 2002.[Eva10] L. C. Evans.
Partial Differential Equations , volume 19 of
Graduate Studies inMathematics . American Mathematical Society, Providence, RI, second edition,2010.[Gal15] D. Gallistl. Stable splitting of polyharmonic operators by generalized Stokessystems.
INS Preprint 1529 , Institut für Numerische Simulation, Germany,2015.[GN11] Th. Gudi and M. Neilan. An interior penalty method for a sixth-order ellipticequation.
IMA J. Numer. Anal. , 31(4):1734–1753, 2011.[Kin89] J. R. King. The isolation oxidation of silicon: the reaction-controlled case.
SIAM J. Appl. Math. , 49(4):1064–1080, 1989.[LM72] J.-L. Lions and E. Magenes.
Non-homogeneous Boundary Value Problems andApplications. Vol. I . Die Grundlehren der mathematischen Wissenschaften,Band 181. Springer-Verlag, New York-Heidelberg, 1972.[Mor68] L.S.D. Morley. The triangular equilibrium element in the solution of platebending problems.
Aeronaut.Quart. , 19:149–169, 1968.[Neč67] J. Nečas.
Les Méthodes Directes en Théorie des Équations Elliptiques . Massonet Cie, Éditeurs, Paris; Academia, Éditeurs, Prague, 1967.[Rud76] W. Rudin.
Principles of Mathematical Analysis . McGraw-Hill Book Co., NewYork-Auckland-Düsseldorf, third edition, 1976.[Sch15] M. Schedensack. A new generalization of the P non-conforming FEM to higherpolynomial degrees. 2015. Preprint, arXiv:1505.02044.[Sch16] M. Schedensack. Mixed finite element methods for linear elasticity and theStokes equations based on the Helmholtz decomposition. ESAIM Math. Model.Numer. Anal. , http://dx.doi.org/10.1051/m2an/2016024, 2016.[Ste08] R. Stevenson. The completion of locally refined simplicial partitions created bybisection.
Math. Comp. , 77(261):227–241, 2008.[SZ90] L. R. Scott and S. Zhang. Finite element interpolation of nonsmooth functionssatisfying boundary conditions.
Math. Comp. , 54(190):483–493, 1990.[Vee14] A. Veeser. Approximating gradients with continuous piecewise polynomial func-tions.
Foundations of Computational Mathematics , pages 1–28, 2014.26Ver96] R. Verfürth.
A Review of a Posteriori Error Estimation and Adaptive Mesh-Refinement Techniques . Advances in numerical mathematics. Wiley, 1996.[WX13] M. Wang and J. Xu. Minimal finite element spaces for m -th-order partialdifferential equations in R n . Math. Comp. , 82(281):25–43, 2013.[Žen70] A. Ženíšek. Interpolation polynomials on the triangle.