FIRST-ORDER SYSTEM LEAST SQUARES ON CURVED
BOUNDARIES: LOWEST-ORDER RAVIART–THOMAS ELEMENTS∗
FLEURIANNE BERTRAND†, STEFFEN M ¨UNZENMAIER‡, AND GERHARD STARKE‡ Abstract. The effect of interpolated edges of curved boundaries on Raviart–Thomas finite element approximations is studied in this paper in the context of first-order system least squares methods. In particular, it is shown that an optimal order of convergence is achieved for lowest-order elements on a polygonal domain. This is illustrated numerically for an elliptic boundary value problem involving circular curves. The computational results also show that a polygonal approximation is not sufficient to achieve convergence of optimal order in the higher-order case.
Key words. interpolated boundaries, Raviart–Thomas spaces, first-order system least squares AMS subject classifications. 65N30, 65N50
DOI. 10.1137/13091720X
1. Introduction. Elliptic boundary value problems arising in science and en-gineering are often posed on domains with curved boundaries. For low-order finite elements it is usually sufficient to use a polygonal approximation of the boundary in order to ensure that the finite element order of convergence is not affected by the inaccurate representation of the boundary. For higher-order finite elements a more accurate representation of the boundary is required which is usually achieved by isoparametric finite elements; i.e., the same finite element space is used once more for the parametrization of a more accurate domain approximation. This is well-known in the case of standard nodal finite element spaces, and a complete theory can be found in many finite element books (see, e.g., [4, 6]).
Finite element methods which simultaneously approximate scalar (potential) and vector (flux) unknowns are also widely used due to the fact that improved approxima-tion of fluxes is sometimes a desirable goal. Mixed and hybrid finite elements of saddle point structure [3] and first-order system least squares [2] are two popular approaches in this direction. Favorable conservation properties are among the advantages of the first of these classes of methods while the latter approach possesses an inherent error control and allows for simplified handling of coupling conditions. On the other hand, properties of either of these approaches are sometimes inherited by the other one due to closeness properties established in [5]. The treatment of curved boundaries in the context of edge- or face-based finite elements where Neumann boundary conditions are imposed on the normal flux is, however, rather seldom mentioned in the liter-ature. A thorough analysis of the effect of polygonal boundary approximation on the Raviart–Thomas finite element approximations seems to be missing. The same appears to be true for the convergence analysis of parametric Raviart–Thomas ele-ments in the higher-order case. General remarks in this direction can be found in the
∗Received by the editors April 16, 2013; accepted for publication (in revised form) January 16,
2014; published electronically April 10, 2014. This work was supported by the German Research Foundation (DFG) under grant STA402/10-1.
http://www.siam.org/journals/sinum/52-2/91720.html
†Institut f¨ur Angewandte Mathematik, Leibniz Universit¨at Hannover, 30167 Hannover, Germany
(bertrand@ifam.uni-hannover.de).
‡Fakult¨at f¨ur Mathematik, Universit¨at Duisburg-Essen, 45127 Essen, Germany (steffen.
muenzenmaier@uni-due.de, gerhard.starke@uni-due.de). 880
standard references on the subject (see, e.g., [3, 10]). The investigation of parametric Raviart–Thomas spaces is deferred to this paper’s sequel.
The purpose of this paper is to provide a convergence analysis of Raviart–Thomas finite elements on curved domains in the context of first-order system least squares formulations. This includes showing optimal order of convergence in the lowest-order case if a polygonal approximation of the boundary is used. It also includes providing numerical evidence of the fact that this is not sufficient in the higher-order case.
For the purpose of exposition, we perform several simplifications. First of all, we restrict ourselves to the two-dimensional case but note that the main results are simi-larly valid in three dimensions. A more severe issue is our restriction to homogeneous boundary conditions. Inhomogeneous boundary conditions for the normal component of the approximated fields cause an additional error which can be controlled in a similar manner using the results derived in this paper. This issue is closely related to the effective treatment of interface conditions in a first-order system least squares approach (see [11]). In this context, the construction of appropriate boundary func-tionals along the lines of [12] for curved boundaries is of interest. A related approach for H1-based first-order system least squares methods is presented in [13].
Throughout this paper, Ω⊂ R2is a bounded domain with a Lipschitz continuous
and piecewise C2 boundary Γ = ∂Ω which is completely covered by a triangulation
Th. Such a triangulation needs to admit curved triangles close to the boundary Γ, in general. Our assumption is that each curved triangle T ∈ Th has the form depicted in Figure 1, i.e., we assume that each element of Th has at most one curved edge. In particular, this means that all points where two smooth boundary segments meet are also vertices of the triangulation. We are concerned with the solution of variational problems in the space
(1.1) HΓ(div, Ω) ={v ∈ H(div, Ω) : n · v, q0,Γ= 0 for all q∈ H1/2(Γ)} . In the definition of HΓ(div, Ω), n denotes the outward pointing normal vector on Γ (which is well-defined due to our smoothness assumptions). The L2(Γ) inner product
·, ·0,Γ is well-defined for traces n· v of functions v ∈ H(div, Ω) if it is regarded as
duality pairing (cf. [8, sect. I.2]).
For the finite element approximation, Ω is approximated by a domain Ωh with polygonal boundary Γh. Associated with Ωh is a triangulation Th which consists completely of straight triangles. For the finite element space Vh, defined with respect to Th, the piecewise polynomial functions vh ∈ Vh may be extended to T ∈ Th by using its polynomial representation on the corresponding element T ∈ Th. However, Vh HΓ(div, Ω), in general, since n· vh does not vanish on Γ for vh∈ Vh. Instead, the condition nh· vh = 0 on Γh may be imposed on the space Vh where nh is the outward normal with respect to Γh. It is necessary to keep track of the effect that this inexactness of the boundary condition has on the overall accuracy of our finite element approximation. For simplicity, our results are formulated for quasiuniform triangulations.
Our focus is on the Raviart–Thomas elements which we define here for a domain Ωh with polygonal boundary Γh. On such a polygonal domain, the Raviart–Thomas space of degree k≥ 0 is given by
(1.2) Vkh= vh∈ HΓh(div, Ωh) : vh|T = qk(x1, x2) rk(x1, x2) + x1 x2 sk(x1, x2) with polynomials qk, rk, and sk of degree k. For the moment, however, we may stick to the case k = 0.
Fig. 1. Polygonal approximation of a curved triangle.
In the next section, an estimate for the normal flux on interpolated boundaries associated with Raviart–Thomas elements is derived. Section 3 shows how this result can be used for the error analysis of a first-order system least squares finite element approximation. The theoretical result is then illustrated by a computational example in section 4. Finally, numerical evidence is provided for the fact that optimal order convergence is not achieved on a polygonal approximate domain in the higher-order case.
2. An estimate for the normal flux on interpolated boundaries. In this section, the subspace V0
h of Raviart–Thomas elements of lowest order on Th with
boundary conditions nh· vh = 0 on Γh is considered. As described above, piecewise polynomial functions vh∈ V0
hmay be extended to Ω. However, Vh0 is not a subspace
of HΓ(div, Ω), in general, for a curved boundary Γ = ∂Ω. The purpose of this section is the derivation of an estimate for the approximated normal flux vh· n on Γ for vh ∈ V0h. To this end, we will use the following abbreviated notation. We simply write a(ξ) b(ξ) if there exists a positive constant C such that a(ξ) ≥ Cb(ξ) holds for all admissible ξ. Further, we write a(ξ) b(ξ) if both a(ξ) b(ξ) and b(ξ) a(ξ) hold. We will repeatedly make use of norm equivalence relations of the form
(2.1) vh20,T hTvh20,∂T, div vh20,T hTdiv vh20,∂T
which hold for all vh ∈ V0h due to a scaling argument and the fact that the restric-tion of the considered space to an element T is finite-dimensional. With the same argument, (2.2) vh20,Ω h vh 2 0,Ω, div vh20,Ωh div vh 2 0,Ω
holds for all vh∈ V0
h.
Theorem 2.1. Assume that Γ = ∂Ω is a piecewise C2 curve. Then, (2.3) |n · vh, q0,Γ| hvh20,Ω+div vh20,Ω1/2q1/2,Γ
holds for all vh∈ V0
h and q∈ H1/2(Γ).
Proof. We consider a fixed curved boundary edge ΓT and the adjacent triangle
T ∈ Th. The coordinate system may be shifted and rotated in such a way that the straight edge ET (cf. Figure 1) is located on the interval [0, hT] of the x1-axis. With respect to these coordinates, ΓT may be parametrized as
(2.4) ΓT = ξ η(ξ) : 0≤ ξ ≤ hT .
The smoothness assumption (C2 parametrization of Γ
T) and the fact that η(0) =
η(hT) = 0 implies that there exists a ¯ξ ∈ [0, hT] such that η( ¯ξ) = 0 and that |η(ξ)| h
T and|η(ξ)| h2T holds for all ξ ∈ [0, hT]. Since nh· vh= 0 on ET,
vh(x) = α 0 + x1 x2 β
holds for all x = (x1, x2)∈ T with α, β∈ R. This implies that on ΓT,
n· vh= 1 (1 + η(ξ)2)1/2 −η(ξ) 1 · α + ξ β η(ξ) β
is satisfied which, for q∈ L2(Γ
T), leads to n · vh, q0,ΓT = ΓT (n· vh) q ds = hT 0 (−η(ξ)(α + ξ β) + η(ξ) β) q(ξ, η(ξ)) dξ hT 0 hT |α + ξ β| + h2T|β||q| dξ .
Using the fact that, on ET, vh(ξ, 0) = α + ξ β 0 and div vh≡ 2β holds, one obtains
n · vh, q0,ΓT hTvh0,ET + h2Tdiv vh0,ETq0,ΓT hTvh20,ET +div vh20,ET
1/2q 0,ΓT ≤ hTvh20,∂T+div vh20,∂T 1/2q 0,ΓT ≤ h1/2 T vh20,T+div vh20,T 1/2 q0,ΓT , (2.5)
where (2.1) is used in the last inequality. Summing over all edges E⊂ Γhleads to
n · vh, q0,Γ h1/2 vh20,Ωh+div vh 2 0,Ωh 1/2q 0,Γ h1/2v h20,Ω+div vh20,Ω 1/2q 0,Γ. (2.6)
For q∈ H1(Γ), integration by parts (note that η(0) = η(h T) = 0) leads to n · vh, q0,ΓT = hT 0 (−η(ξ)(α + ξ β) q(ξ, η(ξ)) + η(ξ) β q(ξ, η(ξ))) dξ = hT 0 η(ξ)(α + ξ β) d dξq(ξ, η(ξ)) + 2η(ξ) β q(ξ, η(ξ)) dξ ≤ hT 0 η(ξ)2(α + ξ β)2dξ 1/2 h T 0 d dξq(ξ, η(ξ)) 2 dξ 1/2 + 2 hT 0 η(ξ)2β2dξ 1/2 h T 0 q(ξ, η(ξ))2dξ 1/2 , which implies
n · vh, q0,ΓT h2T(vh0,ET|q|1,ΓT +div vh0,ETq0,ΓT) h2
T
vh20,ET +div vh20,ET
1/2|q|2 1,ΓT +q20,ΓT 1/2 ≤ h2 T vh20,∂T+div vh20,∂T 1/2q 1,ΓT h3/2 T vh20,T+div vh20,T 1/2 q1,ΓT . (2.7)
Summing over all edges E⊂ Γh,
n · vh, q0,Γ h3/2vh0,Ω2 h+div vh20,Ωh1/2q1,Γ h3/2v h20,Ω+div vh20,Ω 1/2q 1,Γ (2.8)
is obtained. Finally, the result of Theorem 2.1 follows from the fact that H1/2(Γ) is
an interpolation space of type 1/2 for L2(Γ) and H1(Γ) (see [1, Theorem 7.23]).
Remark. Generalizations of (2.3) to triangulations which are not quasiuniform
are of interest in the context of local refinement. In that case, an estimate of the form
|n · vh, q0,Γ| vh20,Ω+div vh20,Ω 1/2 ΓT⊂Γ hTq21/2,Γ T 1/2
would be used instead of (2.3).
3. First-order system least squares with curved boundaries. The result of Theorem 2.1 is now applied to the first-order system least squares formulation using Raviart–Thomas elements for the flux approximation. For simplicity, we restrict our analysis to the first-order system corresponding to Poisson’s equation
div u = f, u +∇p = 0 (3.1)
with given f ∈ L2(Ω) subject to Neumann boundary conditions n· u = 0 on Γ = ∂Ω.
For the finite element approximation of the flux, the previously described lowest-order Raviart–Thomas elements V0
h are used. The pressure p is approximated by
conforming piecewise linear elements on a triangulation of Ωh. The least squares functional associated with this problem is
(3.2) F(u, p) = div u − f20,Ω+u + ∇p20,Ω
to be minimized with respect to u∈ HΓ(div, Ω) and p∈ H1(Ω), and for the approxi-mation on Ωh we use (3.3) Fh(u, p) =div u − fh20,Ω h+u + ∇p 2 0,Ωh.
In order to have a unique minimizer of (3.2) and (3.3) an additional condition on p is required, e.g., (p, 1)0,Ω= 0. With the definition of the subspace
(3.4) H˙1(Ω) ={q ∈ H1(Ω) : (q, 1)0,Ω= 0} , the Poincar´e–Friedrichs inequality
(3.5) q20,Ω≤ CF∇q20,Ω for all q∈ ˙H1(Ω)
holds (cf. [4, sect. II.1]). The finite element space for the pressure approximation is then given by
(3.6) Q1h={qh∈ ˙H1(Ω) : qh|T is a polynomial of degree 1 for each T ∈ Th} . The trace inequality
(3.7) q21/2,Γ≤ CTq21,Ω for all q∈ H1(Ω)
(cf. [8, Theorem 1.5]) will also be used in what follows. More precisely, combined with (3.5), (3.7) implies
(3.8) q21/2,Γ≤ CT(1 + CF)∇q21,Ω=: CT∇q21,Ω for all q∈ ˙H1(Ω) , which is the form of the trace inequality used in the proof of Theorem 3.3.
In practice, the constraint (qh, 1)0,Ω = 0 will be replaced in Q1
h by the
imple-mentable condition (qh, 1)0,Ωh = 0. This difference, however, has no effect on our error estimates since the two solutions only differ by a constant on the order of h which can be seen as follows: Let phbe the piecewise linear conforming finite element approximation satisfying ( ph, 1)0,Ωh = 0, then, in order to project phinto Q1
h, one only
needs to subtract the constant ( ph, 1)0,Ω/(1, 1)0,Ω. This constant can be estimated as ( ph, 1)0,Ω= ( ph, 1)0,Ω− ( ph, 1)0,Ωh = ( ph, 1)0,Ω\Ωh− ( ph, 1)0,Ωh\Ω = T ∈ T◦ h ( ph, 1)0, T \T − ( ph, 1)0,T \ T , (3.9)
whereTh◦ denotes those triangles ofThadjacent to the boundary ∂Ω. With an affine functionχ : T ∪ T → R2such that divχ = 1, n · χ = 0 on E
T, and (n× χ, 1)0,ET = 0, using the divergence theorem, we obtain
(3.10) ( ph, 1)0, T \T ≤ ph, n· χ0,∂ T ∩∂Ω +(∇ ph,χ)0, T \T .
The first term on the right-hand side in (3.10) can be estimated as in (2.5) in the proof of Theorem 2.1 leading to
ph, n· χ0,∂ T ∩∂Ω h1/2T χ2 0, T+div χ 2 0, T 1/2 ph0,∂ T ∩∂Ω h1/2 T χ0, T+ hT ph1, T,
where (3.8) was used. Combining this with a reference element calculation which shows thatχ0, T h2
T holds and inserting this into (3.10) leads to
(3.11) ( ph, 1)0, T \T
χ0, T + h3/2T
ph1, T h3/2T ph1, T .
Combined with a similar estimate for ( ph, 1)0,T \ T, we end up with
|( ph, 1)0,Ω| ≤ T ∈ T◦ h |( ph, 1)0, T \T| + |( ph, 1)0,T \ T| T ∈ T◦ h h3/2T ph1, T ≤ h ⎛ ⎝ T ∈ T◦ h hT ⎞ ⎠ 1/2 ph21, T 1/2 h ph1,Ω, (3.12)
where the fact was used that
T ∈ T◦
h
hT =|∂Ω| ,
the length of the boundary.
A connection between the least squares functionalsF and Fh, defined on Ω and Ωh, respectively, is established by an auxiliary mapping Φh. We summarize the main properties of this mapping in the following lemma (cf. [6, sect. 10.4] or [9]).
Lemma 3.1. There exists a mapping Φh : Ω → Ωh such that Φh( T ) = T is
satisfied for all T ∈ Th, T ∈ Th. Moreover, if JΦh denotes the Jacobian of Φh(x),
then JΦh ∈ W∞1(Ω), JΦ−1 h ∈ W 1 ∞(Ω), and (3.13) Φh− idL ∞( T ) h2T , JΦh− IL∞( T ) hT, andJΦ−1h − IL∞( T ) hT
holds for all T∈ Th. For each v : Ω→ R2, ifv : Ω◦
h→ R2 is defined by
(3.14) v(Φ◦ h(x)) = 1
det JΦh(x)JΦh(x)v(x) ,
thenv◦ ∈ HΓh(div, Ωh) if and only if v∈ HΓ(div, Ω), and (3.15) (div◦ v)(Φ◦ h(x)) = 1
det JΦh(x)(div v)(x) for all x∈ Ω
holds (where div denotes differentiation with respect to the new coordinates◦ x =◦ Φh(x)∈ Ωh). Moreover, v◦ ∈ H1(Ωh)2 if and only if v∈ H1(Ω)2.
Proof. For the existence of a mapping Φh satisfying (3.13), see [6, sect. 10.4] and [9]. For better distinction, we denote differentiation with respect tox = Φ◦ h(x) by∇,◦
◦
div, etc. The formula (3.15) is a well-known result and follows from the fact that, for any φ∈ H1 0(Ω) (and corresponding ◦ φ∈ H1 0(Ωh) defined by ◦ φ(Φh(x)) = φ(x)), −(div◦ v,◦ φ)◦ 0,Ωh = (v,◦ ∇◦φ)◦ 0,Ωh = Ωh ◦ v·∇◦φ d◦ x◦ = Ω 1 det JΦh(JΦhv)· (J −T Φh ∇φ) det JΦhdx = Ω v· ∇φ dx = (v, ∇φ)Ω=−(div v, φ)Ω
holds which leads to (3.16) Ω (div◦ v)(Φ◦ h(x)) φ(x) det JΦh(x) dx = Ω (div v)(x) φ(x) dx and implies (3.15). Similarly,
n◦·v,◦ φ◦0,∂Ωh = (div◦ v,◦ φ)◦ 0,Ωh+ (v,◦ ∇◦φ)◦ 0,Ωh
= (div v, φ)0,Ω+ (v,∇φ)0,Ω =n · v, φ0,∂Ω,
which implies thatn◦·v = 0 on ∂Ω◦ h is equivalent to n· v = 0 on ∂Ω. Finally, the chain rule applied to (3.14) gives
(∇◦v)J◦ Φh = 1
det JΦh ([JΦh∂1v , JΦh∂2v] + [(∂1JΦh)v , (∂2JΦh)v])
− 1
(det JΦh)2(JΦhv) [(∂1JΦh) : CofJΦh, (∂2JΦh) : CofJΦh] , which implies the equivalence ofv◦ ∈ H1(Ωh)2 and v∈ H1(Ω)2.
In order to get the optimal order of convergence stated in Theorem 3.3, fh in the discrete functional (3.3) has to be a sufficiently good approximation to f on Ω, in general. Practically, this can be achieved by using the L2(Ω
h)-orthogonal projection
onto piecewise constant functions which is the statement of the following lemma. Lemma 3.2. For Ω ⊂ R2 with curved boundary segments and Ωh ⊂ R2 with
polygonal boundary satisfying the assumptions in the introduction, letPh: L2(Ω
h)→
Zh be the L2(Ω
h)-orthogonal projection onto the space Zh of piecewise constants on
the triangulation Th of Ωh. Then,
(3.17) f − Phf0,Ω hfW1
∞(Ω∪Ωh)
holds for any f ∈ W1
∞(Ω∪ Ωh)
Proof. For each T ∈ Th, orthogonality with respect to L2(T ) implies
Phf|T = (f, 1)0,T (1, 1)0,T and f − Phf 2 0,T =f20,T− (f, 1)2 0,T (1, 1)0,T . On the other hand, for each curved triangle T ∈ Th,
f − Phf2 0, T =f20, T− 2(f, 1)0, T (f, 1)0,T (1, 1)0,T + (1, 1)0, T (f, 1)20,T (1, 1)2 0,T =f2 0, T− (f, 1)2 0, T (1, 1)0, T + (1, 1)0, T (f, 1)0, T (1, 1)0, T − (f, 1)0,T (1, 1)0,T 2 =f − Phf2 0, T + (1, 1)0, T (f, 1)0, T (1, 1)0, T − (f, 1)0,T (1, 1)0,T 2 , (3.18) where ˜Ph: L2(Ω)→ Z
hdenotes the L2(Ω)-orthogonal projection.
For the first term in (3.18) we have, due to the best approximation property,
(3.19) f − Phf2
0, T h 2
Tf21, T .
The second term in (3.18) needs more careful consideration. To this end, we first observe (using, for example, the same parametrization as in the beginning of the proof of Theorem 2.1) that
(3.20) (1, 1)0,T− (1, 1)0, T = |T | − | T| h3T
holds. Second, it is convenient to use the map Φh : Ω→ Ωh introduced in Lemma 3.1. This leads us to (f, 1)0,T − (f, 1)0, T = T (f (Φh(x)) det JΦh(x)− f(x)) dx = T (f (Φh(x))− f(x)) det JΦh(x) dx + T f (x) (det JΦh(x)− 1) dx , (3.21)
which we may further use to get (3.22) (f, 1)0,T− (f, 1)0, T
T|f(Φh
(x))− f(x)| dx + f0, T det JΦh− 10, T.
From the fact that, for all x∈ T , |f(Φh(x))− f(x)| = 1 0 ∇f(x + s(Φh(x)− x)) · (Φh(x)− x) ds ≤ 1 0 |∇f(x + s(Φh(x)− x))| |Φh(x)− x| ds holds, and using (3.13), (3.22) turns into
(3.23) (f, 1)0,T− (f, 1)0, T h3TfW1
∞(T ∪ T )+ h
2
Tf0, T.
Combining (3.20) and (3.23) leads to (f, 1)0, T (1, 1)0, T − (f, 1)0,T (1, 1)0,T = (f, 1)0, T((1, 1)0,T− (1, 1)0, T) (1, 1)0, T(1, 1)0,T + (f, 1)0, T − (f, 1)0,T (1, 1)0,T h−1 T |(f, 1)0, T| + hTfW∞1(T ∪ T )+f0, T hTfW∞1(T ∪ T )+f0, T. (3.24)
Inserting (3.24) and (3.19) into (3.18), (3.25) f − Phf2 0, T h 2 T f2 1, T+ h 2 Tf2W1 ∞(Ω∪Ωh) is obtained. Summing over all triangles finally leads to (3.17).
Theorem 3.3. Assume that the exact solution of the system (3.1) (u, p) ∈
HΓ(div, Ω)× ˙H1(Ω) satisfies p∈ H2(Ω) (implying u∈ H1(Ω)2) and div u∈ H1(Ω).
Moreover, assume that in (3.3) fh =Phf is set, where Ph is the L2(Ω
h)-orthogonal
projection onto the space of piecewise constant functions. Then, if (uh, ph)∈ V0
h×Q1h
denotes its finite element approximation,
(3.26) u − uhdiv,Ω+p − ph1,Ω hu1,Ω+fW1
∞(Ω∪Ωh)+p2,Ω
holds.
Proof. Using the fact that (u, p)∈ HΓ(div, Ω)× ˙H1(Ω) is the exact solution of
the first-order system (3.1), one obtains
F(uh, ph) =div uh− f0,Ω2 +uh+∇ph20,Ω
=div(uh− u)20,Ω+uh− u + ∇(ph− p)20,Ω =div(uh− u)20,Ω+uh− u20,Ω+∇(ph− p)20,Ω
+ 2(1− γ)(uh− u, ∇(ph− p))0,Ω
+ 2γ (n · uh, ph− p0,Γ− (div(uh− u), ph− p)0,Ω) with γ∈ [0, 1] still free to be chosen appropriately. This leads to
F(uh, ph)≥ div(uh− u)20,Ω+ γuh− u20,Ω+ γ∇(ph− p)20,Ω
−γ δdiv(uh− u) 2 0,Ω− γδph− p20,Ω− 2γ |n · uh, ph− p0,Γ| ≥1−γ δ div(uh− u)20,Ω+ γuh− u20,Ω + γ(1− δCF)∇(ph− p)20,Ω− 2γ |n · uh, ph− p0,Γ| , (3.27)
where δ > 0 is still free to be chosen appropriately and the Poincar´e–Friedrichs in-equality (3.5) is used. If this is combined with the inin-equality
2|n · uh, ph− p0,Γ| ≤ 2C1huh20,Ω+div uh20,Ω1/2ph− p1/2,Γ ≤1 εC2h 2u h20,Ω+div uh20,Ω + εC3ph− p21/2,Γ ≤1 εC2h 2u h20,Ω+div uh20,Ω + εC3CT∇(ph− p)20,Ω,
where ε > 0 is still free to be chosen appropriately and where Theorem 2.1 and the trace inequality (3.8) are used, then
F(uh, ph) +γεC2h2 uh20,Ω+div uh20,Ω ≥1−γ δ div(uh− u)20,Ω+ γuh− u20,Ω + γ(1− δCF + εC3CT)∇(ph− p)20,Ω (3.28)
follows. Choosing δ = 1/(3CF), ε = 1/(3C3CT), and γ = δ/2 implies
F(uh, ph)+ h2uh20,Ω+div uh20,Ω
div(uh− u)20,Ω+uh− u20,Ω+∇(ph− p)20,Ω.
(3.29)
Using
uh20,Ω≤ 2u20,Ω+ 2u − uh20,Ω, div uh20,Ω≤ 2div u20,Ω+ 2div(u − uh)20,Ω,
(3.29) implies that
F(uh, ph)+ h2u20,Ω+div u20,Ω
div(uh− u)20,Ω+uh− u20,Ω+∇(ph− p)20,Ω
(3.30)
holds at least for h ≤ h0 with sufficiently small h0 (depending on the constants in (3.28)). Since the left-hand side in (3.30) is bounded from below by a positive constant for h≥ h0, this estimate actually holds for all h.
An upper bound for the least squares functional is obtained based on
F(uh, ph) = T ∈ Th div uh− f20, T +uh+∇ph20, T f − fh20,Ω+ T ∈Th div uh− fh20,T+uh+∇ph20,T =f − fh20,Ω+Fh(uh, ph) h2f2W1 ∞(Ω∪Ωh)+Fh(uh, ph) , (3.31)
which holds due to the fact that the spaces V0
h and Q1h restricted to an individual
element are finite-dimensional and due to Lemma 3.2.
Using the mapping Φh: Ω→ Ωhfrom Lemma 3.1, we define, with respect to the exact solution (u, p)∈ HΓ(div, Ω)× ˙H1(Ω) of (3.1),
◦ u∈ HΓh(div, Ωh) byu(Φ◦ h(x)) = 1 det JΦh(x)JΦh(x)u(x) , ◦ p∈ H1(Ωh) byp(Φ◦ h(x)) = p(x) . (3.32)
Due to the fact that (uh, ph) minimizesFh(uh, ph) in the finite element space Vh×Qh,
Fh(uh, ph)≤ Fh(Rhu,◦ Ihp)◦ =div(◦ Rhu)◦ − fh20,Ω h+Rh ◦ u +∇(I◦ hp)◦ 20,Ω h =Ph(div◦ u)◦ − fh20,Ω h+Rh ◦ u−u +◦ u +◦ ∇◦p◦−∇◦p +◦ ∇(I◦ hp)◦ 20,Ω h =Rhu◦−u◦20,Ωh+u +◦ ∇◦p◦20,Ωh+∇(I◦ hp◦−p)◦ 20,Ωh (3.33)
holds, whereRh : H(div, Ωh)→ V0h and Ih : H1(Ωh)→ Q1h are the corresponding interpolation operators. The approximation results for the Raviart–Thomas spaces (cf. [3, sect. 2.5]) and standard conforming elements (cf. [4, sect. II.6] and [6, sect. 4.4]), respectively, lead to
Rhu◦ −u◦0,Ωh h u◦1,Ωh u1,Ω, ∇(I◦ hp− p)0,Ωh hp◦2,Ωh p2,Ω,
(3.34)
where the rightmost inequalities are due to Lemma 3.1 and to the direct use of the chain rule, respectively. Using the mapping Φhwe also get
u +◦ ∇◦p◦20,Ω h = Ωh ◦ u +∇◦p◦ 2 dx◦ = Ω 1 det JΦhJΦhu + J −T Φh ∇p 2 det JΦhdx = Ω 1 det JΦhJΦh− I u +JΦ−T h − I ∇p 2 det JΦhdx,
which implies u +◦ ∇p◦20,Ωh 1 det JΦhJΦh− I 2 L∞(Ω) u2 0,Ω+JΦ−1h − IL∞(Ω)∇p20,Ω h2u2 0,Ω+∇p20,Ω , (3.35)
where we used (3.13) again. Inserting (3.34) and (3.35) into (3.33) leads to (3.36) Fh(uh, ph) h2u21,Ω+p22,Ω.
The result finally follows from combining (3.30), (3.31), (3.36), and (3.5).
Remark. For polygonal domains, Theorem 3.3 reduces to the convergence result
obtained in [7].
Remark. We restrict ourselves to homogeneous boundary conditions n· u = 0
throughout this paper. The more general case n· u = g with g ∈ H−1/2(Γ) could be treated in the following way: Write u = uN + u0 with n· uN = g on Γ and
uh= uNh + u0
hwith u0h∈ V0h (i.e., n· u0h= 0 on Γh). Then, the “homogeneous” error
term u0− u0
h can be treated by the techniques developed above while the remainder
uN− uNh is a question of approximability of g◦ Φ−1h by piecewise constant functions on Γh.
4. Computational results. This section starts with an illustration of our con-vergence result in Theorem 3.3 by numerical computations. The numerical tests are performed for the first order system (3.1) on the unit disk Ω ={(x1, x2) : x2
1+ x22< 1}
with boundary conditions n· u = 0 on ∂Ω. For the first example, the right-hand side
f is chosen in such a way that the exact solution of (3.1) is given by
(4.1) u(x1, x2) = 4x1(1− x2 1− x22) 4x2(1− x2 1− x22) , p(x1, x2) = (1− x21− x22)2−1 3. This leads to f (x1, x2) = 8− 16(x21+ x22), the exact solution is shown in Figure 2.
−1 −0.5 0 0.5 1 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 ||u|| 2 0 0.5 1 1.5 (a)∇p = −u −1 −0.5 0 0.5 1 −1 −0.5 0 0.5 1 −0.2 0 0.2 0.4 0.6 (b)p
Fig. 2. Exact solution (4.1).
The computed values of the least squares functionalFh(uh, ph) are shown on the left in Figure 3. For the V0
h/Q1h combination a convergence rate appears which is
inversely proportional to the number of unknowns N . The same convergence rate
103 104 105 10−8 10−6 10−4 10−2 100 Degrees of freedoms
Least Squares Functional
1 1 1.8 1 RT0−P1 RT1−P2 103 104 105 10−8 10−7 10−6 10−5 10−4 10−3 10−2 10−1 Degrees of freedoms ||u−u h || 2 + || ∇ p− ∇ ph || 2 1 1 1.8 1 RT0−P1 RT1−P2
Fig. 3. Convergence rates for the first-order system with exact solution (4.1).
is observed for the squares of the error normsu − uh2
0,Ωh and ∇(p − ph)20,Ωh in the right graph of Figure 3. Due to the fact that div(u − uh)20,Ω
h is part of the functional Fh(uh, ph) and due to (3.5), this implies the same order of convergence (inversely proportional to N ) foru − uh2div,Ω
h andp − ph
2
1,Ωh. Since the number of degrees of freedom is proportional to h−2 in two dimensions, this confirms the predictions by Theorem 3.3, which is the optimal order of convergence achievable. Figure 3 also shows that the convergence rate for the V1
h/Q2h combination is quite a bit faster. 103 104 105 106 10−8 10−7 10−6 10−5 10−4 10−3 10−2 Degrees of freedoms
Least Squares Functional
1 1 1.5 1 RT0−P1 RT1−P2
Fig. 4. Convergence rates for the first order system (3.1) withf(x1, x2) = sin(x1). However, this example does not provide a clear indication of whether the optimal
order of approximation achievable is actually attained or not in the V1
h/Q2h case.
Therefore, we consider a second example, again on the unit disk using the same boundary conditions and normalizing constraint as above. Here, the right-hand side is given by f (x1, x2) = sin(x1). Figure 4 clearly shows that the order of convergence is no longer optimal for the V1h/Q2hcombination. A closer look at the proof of Theorem 3.3 shows where modifications are needed in order to obtain the optimal convergence behavior
(4.2) u − uh2div,Ω+p − ph21,Ω h4div u22,Ω+u22,Ω+p23,Ω.
For the required improvement of (3.30), a higher-order version of Theorem 2.1 is needed. In order to improve (3.36) to the higher-order case, the interpolation opera-torsRhandIhneed to be adjusted. These modifications will be presented in a future paper based on suitable parametric finite element spaces.
We close our presentation of numerical experiments with a three-dimensional version of (4.1) in the unit ball Ω ={(x1, x2, x3) : x12+ x22+ x23 < 1} (see Figure 5).
The results in Figure 6 show that the functional as well as the square norm of the error behave in an optimal way proportionally to h2∼ N−2/3in the lowest-order case.
In the higher-order case (RT1combined with P2 elements) the convergence behavior is similar to the two-dimensional situation as well.
Fig. 5. Exact solution and triangulation on the unit ball.
5. Conclusions. Optimal order convergence of a first-order system least squares method using lowest-order Raviart–Thomas combined with linear conforming ele-ments is shown for domains with curved boundaries. In particular, an estimate for
103 104 105
102
103
104
Degrees of freedoms
Least Squares Functional
1.3 1 2/3 1 RT0−P1 RT1−P2 103 104 105 102 103 104 Degrees of freedoms ||u−u h || 2 + || ∇ p− ∇ ph || 2 1 1 1.3 2/3 RT0−P1 RT1−P2
Fig. 6. Convergence rates for the first-order system in the three-dimensional unit ball.
the normal flux associated with Raviart–Thomas elements on interpolated bound-aries is derived. Moreover, the convergence proof is based on mapping the polygonal approximate domain onto the original domain. Our computational results confirm the theory and suggest that a polygonal domain approximation is not sufficient in the higher-order case. An appropriate approach based on parametric finite element spaces will be presented in a future paper.
Acknowledgment. We would like to thank the anonymous referees for their careful reading and valuable suggestions.
REFERENCES
[1] R. A. Adams and J. F. Fournier, Sobolev Spaces, 2nd ed., Elsevier/Academic Press, Amster-dam, 2003.
[2] P. B. Bochev and M. D. Gunzburger, Least-Squares Finite Element Methods, Springer, New York, 2009.
[3] D. Boffi, F. Brezzi, and M. Fortin, Mixed Finite Element Methods and Applications, Springer, Heidelberg, 2013.
[4] D. Braess, Finite Elements: Theory, Fast Solvers, and Applications in Solid Mechanics, 2nd ed., Cambridge University Press, Cambridge, 2001.
[5] J. Brandts, Y. Chen, and J. Yang, A note on least-squares mixed finite elements in relation to standard and mixed finite elements, IMA J. Numer. Anal., 26 (2006), pp. 779–789. [6] S. C. Brenner and L. R. Scott, The Mathematical Theory of Finite Element Methods,
3rd ed., Springer, New York, 2008.
[7] Z. Cai, R. Lazarov, T. A. Manteuffel, and S. F. McCormick, First-order system least squares for second-order partial differential equations: Part I, SIAM J. Numer. Anal., 31 (1994), pp. 1785–1799.
[8] V. Girault and P.-A. Raviart, Finite Element Methods for Navier-Stokes Equations, Springer-Verlag, Berlin, 1986.
[9] M. Lenoir, Optimal isoparametric finite elements and error estimates for domains involving curved boundaries, SIAM J. Numer. Anal., 23 (1986), pp. 562–580.
[10] P. Monk, Finite Element Methods for Maxwell’s Equations, Oxford University Press, New York, 2003.
[11] S. M¨unzenmaier and G. Starke, First-order system least squares for coupled Stokes–Darcy flow, SIAM J. Numer. Anal., 49 (2011), pp. 387–404.
[12] G. Starke, Multilevel boundary functionals for least-squares mixed finite element methods, SIAM J. Numer. Anal., 36 (1999), pp. 1065–1077.
[13] R. Stevenson, A note on first order system least squares with inhomogeneous boundary con-ditions, IMA J. Numer. Anal., to appear.
