# Quasi-Optimality of Adaptive Mixed FEMs for Non-selfadjoint Indefinite Second-Order Linear Elliptic Problems

Carsten Carstensen, Asha K. Dond and Hella Rabus

# Abstract

The well-posedness and the a priori and a posteriori error analysis of the lowest-order Raviart–Thomas mixed finite element method (MFEM) has been established for non-selfadjoint indefinite second-order linear elliptic problems recently in an article by Carstensen, Dond, Nataraj and Pani (Numer. Math., 2016). The associated adaptive mesh-refinement strategy faces the difficulty of the flux error control in H(div,Ω) and so involves a data-approximation error f-Π0f in the L2 norm of the right-hand side f and its piecewise constant approximation Π0f. The separate marking strategy has recently been suggested with a split of a Dörfler marking for the remaining error estimator and an optimal data approximation strategy for the appropriate treatment of f-Π0fL2(Ω). The resulting strategy presented in this paper utilizes the abstract algorithm and convergence analysis of Carstensen and Rabus (SINUM, 2017) and generalizes it to general second-order elliptic linear PDEs. The argument for the treatment of the piecewise constant displacement approximation uRT is its supercloseness to the piecewise constant approximation Π0u of the exact displacement u. The overall convergence analysis then indeed follows the axioms of adaptivity for separate marking. Some results on mixed and nonconforming finite element approximations on the multiply connected polygonal 2D Lipschitz domain are of general interest.

## 1 Introduction

The convergence analysis of adaptive mixed finite element methods (AMFEM) stated in [6, 7, 11, 8] for the Laplacian is completed in this paper for non-selfadjoint indefinite second-order linear elliptic problems via separate marking with the axioms from [8]. Given a right-hand side fL2(Ω) and piecewise constant coefficients in a (possibly multiply connected) bounded polygonal Lipschitz domain Ω2, the general second-order linear elliptic PDE seeks uH01(Ω) such that

(1.1)u:=-div(𝐀u+u𝐛)+γu=fin Ω.

The coefficients 𝐀, 𝐛, and γ are all piecewise constant functions and the symmetric matrix 𝐀 is assumed to be positive definite with universal positive eigenvalue bounds from below and above. The entire paper solely assumes that :H01(Ω)H-1(Ω) is injective (i.e., has a trivial kernel); then it is bijective with a bounded inverse and satisfies a global inf-sup condition. The flux variable 𝐩=-(𝐀u+u𝐛) with 𝐛=𝐀-1𝐛 allows to recast problem (1.1) as a first-order system: Seek uH01(Ω) such that

(1.2)𝐀-1𝐩+u𝐛+u=0anddiv𝐩+γu=fin Ω.

The mixed formulation of (1.1) seeks (𝐩,u)H(div,Ω)×L2(Ω) such that

(1.3)

(𝐀-1𝐩+u𝐛,𝐪)L2(Ω)-(div𝐪,u)L2(Ω)=0,
(div𝐩,v)L2(Ω)+(γu,v)L2(Ω)=(f,v)L2(Ω)

for all 𝐪H(div,Ω) and for all vL2(Ω). The well-posedness of (1.3) has been studied in [4, Theorem 2.1] using the equivalence of the weak formulation of (1.1) and the mixed formulation (1.3). The mixed finite element discretization of (1.3) utilizes the piecewise constant functions P0(𝒯) on 𝒯 and the lowest-order Raviart–Thomas finite element space RT0(𝒯)H(div,Ω) and seeks (𝐩RT,uRT)RT0(𝒯)×P0(𝒯) such that

(1.4)

(𝐀-1𝐩RT+uRT𝐛,𝐪RT)L2(Ω)-(div𝐪RT,uRT)L2(Ω)=0,
(div𝐩RT,vRT)L2(Ω)+(γuRT,vRT)L2(Ω)=(Π0f,vRT)L2(Ω)

for all 𝐪RTRT0(𝒯) and for all vRTP0(𝒯). The well-posedness of (1.4) has been established in [4, Theorem 4.2] under the additional assumption that the initial mesh-size is sufficiently small.

The convergence and quasi-optimality of adaptive finite element methods for linear symmetric elliptic problems has been discussed in the literature [1, 2, 5, 7, 8, 9, 17, 18, 20, 21] and the references mentioned therein. For the non-symmetric case 𝐛0 and for adaptive conforming FEMs, the optimal convergence rates have been established in [15] with a collective marking strategy based on a posteriori error estimation, and – in contrast to the results in [16] – without the interior node property for the refinement. The convergence of adaptive nonconforming FEMs for the non-symmetric and indefinite problem has been derived in [10] with a different separate marking strategy. The adaptive mixed FEM (1.4) has been considered in [13] with a combined norm of the L2 norm in the flux error and the L2 norm in the displacement error. Their quasi-optimality analysis is based on some nonstandard adaptive separate marking scheme and a special relation between the mixed FEM and nonconforming schemes.

This paper develops the quasi-optimality of adaptive MFEMs with the natural H(div) norm, that is, the combination of a flux error in the H(div) norm and the displacement error in the L2 norm via the axioms for separate marking from [8]. The total adaptive estimator is the sum of the residual-based error estimator η(𝒯) and the data approximation error μ(𝒯). Given a parameter κ>0, the separate marking scheme runs either the Dörfler marking [14] on η(𝒯) if μ2(𝒯)κη2(𝒯) or an optimal data approximation algorithm as in [7, 17, 18] to reduce μ(𝒯). The main challenge is the proof of the axioms of discrete reliability and quasi-orthogonality for the non-symmetric mixed problem.

The first intermediate solution concerns the discrete flux approximation with a prescribed divergence on the coarse triangulation in the finer Raviart–Thomas space and a generalization of the corresponding design from [8, 7, 11]. The second intermediate solution is the integral mean Π0u of the exact displacement u and its supercloseness

(1.5)Π0u-uRTChmaxα(𝐩-𝐩RTH(div,Ω)+u-uRT)

in the proof of quasi-orthogonality (A4ϵ) below. In fact, this difficulty does not arise in [8] for 𝐛0. At first glance, the extended Marini-type identity [4, equation (4.3)] states that uCR is close to Π0uCR for some associated Crouzeix–Raviart solution uCR, which is superclose to u by L2-duality arguments with α>0 from reduced elliptic regularity. This argument, however, leads to (1.5) up to an additional term h𝒯div𝐩RT, which is not included in the error estimator η utilized in this paper. In fact, η solely involves the term h𝒯2div𝐩RT with a higher power of the localized mesh-size h𝒯.

The remaining parts of the paper are organized as follows. Section 2 establishes the notation, the a posteriori error estimators, the adaptive algorithm with separate marking (Safem), and recalls the axioms of adaptivity (A1)(A4), (B1)(B2), and quasimonotonicity (QM) with the optimal convergence rates from [8]. Section 3 starts with the proof of stability (A1) and reduction (A2) for the error estimators and distance functions at hand. Section 4 is devoted to the discrete reliability based on a discrete Helmholtz decomposition in 2D for multiply connected domains. Section 5 verifies the quasi-orthogonality based on (1.5) with a direct proof in Lemma 1. Numerical experiments in Section 6 investigate the condition on sufficiently small parameters such as the bulk parameter and the mesh-size for optimal convergence rates.

The presentation is laid out for two-dimensional polygonal domains and the lowest-order case. The coefficients are assumed piecewise constant for simplicity to avoid extra perturbation terms as in [9]. The generalization to 3D may follow the lines of this paper and replaces the discrete Helmholtz decomposition as in [7] for a simply connected domain. The analysis of higher-order finite element approximations requires a new argument for the stability of the discrete problems in that [4] and part of the proofs in this paper utilize the equivalence to the Crouzeix–Raviart nonconforming FEM, which is open in 3D for higher polynomial degrees.

Standard notation on Lebesgue and Sobolev spaces such as L2(Ω), H01(Ω), and H(div,Ω) and their norms with :=L2(Ω) and :=L(Ω) apply throughout the paper. The notation AB abbreviates ACB for a mesh-size independent generic constant C>0, which may depend on the domain Ω and the shape but not the size of the triangles of the corresponding shape-regular triangulation. The constant C may also depend on the coefficients 𝐀, 𝐛, and γ through lower and upper bounds of the positive eigenvalues of 𝐀 and the upper bound 𝐛+γ for the remaining coefficients. Furthermore, AB abbreviates AB and BA. The context depending symbol || denotes the area of a domain, the length of an edge, the counting measure (cardinality) of a set, the absolute value of a real number, or the Euclidean length of a vector.

## 2 Preliminaries

This section first introduces the necessary notation for the definition and analysis of adaptive algorithms with separate marking. The axioms of adaptivity from [8] are slightly simplified to match the setting of this paper.

### 2.1 Notation

Let 𝒯 be an admissible triangulation of the bounded polygonal domain Ω and let 𝕋(𝒯) be the set of all admissible triangulations refined from 𝒯 by newest vertex bisection [20]. Let (T) denote the set of the three edges of a triangle T𝒯, let (resp. (Ω) and (Ω)) denote the set of all (resp. interior and boundary) edges in the triangulation 𝒯 and let 𝒩 be the set of its vertices.

Let hmax:=maxT𝒯hT denote the maximal local mesh-size hT:=|T|12 and let νE and τE are the unit normal and tangential vectors along E(T) of T𝒯. The jump [𝐪]E:=𝐪|T+-𝐪|T- of 𝐪 is defined across an interior edge E shared by the two triangles T+ and T-, which form the edge patch ωE¯. For any boundary edge EΩ, let ωE=T+ denote the interior of the triangle T+=ωE¯ with the edge E(ωE¯) and the jump []E reduces to the trace (that is, the exterior jump contribution vanishes according to the homogeneous boundary condition). For vH1(Ω;) and 𝚽:=(ϕ1,ϕ2)H1(Ω;2), the curl and gradient operators read

Curlv=(-vy,vx),curlΦ=ϕ1y-ϕ2x,v=(vx,vy).

The piecewise gradient NC acts as (NCv)|T=(v|T) for all T𝒯. Let Pr denote the algebraic polynomials of degree at most r and set

Pr(𝒯):={vL2(Ω):v|TPr(T) for all T𝒯},S1(𝒯):=P1(𝒯)C(Ω).

Let Π0 denote the piecewise L2 projection onto P0(𝒯) with respect to the shape-regular triangulation 𝒯. The associated nonconforming Crouzeix–Raviart and lowest-order Raviart–Thomas mixed finite element spaces read

CR1(𝒯):={vP1(𝒯):v is continuous in all midpoints mid(E) of edges E},
CR01(𝒯):={vCR1(𝒯):v(mid(E))=0 for all E(Ω)},
RT0(𝒯):={𝐪H(div,Ω):T𝒯𝐜2d𝐱T,𝐪(𝐱)=𝐜+d𝐱}.

### 2.2 A Posteriori Error Control

The a posteriori error control for problem (1.4) has been established in [4, Theorem 5.5]: Given the unique solutions (𝐩,u) to (1.3) and (𝐩RT,uRT) to (1.4), for small initial mesh-size h0, the equivalence

(2.1)𝐩-𝐩RTH(div,Ω)2+u-uRT2σ2(𝒯):=η2(𝒯)+μ2(𝒯)

holds for the (squared) error estimator η(𝒯) and data approximation μ(𝒯) defined by

(2.2)η2(𝒯):=T𝒯η2(𝒯,T)with
η2(𝒯,T):=|T|12E(T)[𝐀-1𝐩RT+uRT𝐛]EτEL2(E)2+|T|𝐀-1𝐩RT+uRT𝐛L2(T)2

and

(2.3)μ2(𝒯):=T𝒯μ2(T)with μ2(T):=f-Π0fL2(T)2.

Given the discrete solution (𝐩RT,uRT) and (𝐩RT^,uRT^) with respect to the admissible triangulation 𝒯 and its refinement 𝒯^𝕋(𝒯), respectively, the distance function reads

(2.4)δ2(𝒯,𝒯^):=𝐩RT^-𝐩RTH(div,Ω)2+uRT^-uRT2

with the weighted norm from

𝐩RT^-𝐩RTH(div,Ω)2:=𝐀-12(𝐩RT^-𝐩RT)2+div(𝐩RT^-𝐩RT)2.

### 2.3 Safem – The Adaptive Algorithm with Separate Marking

The separate marking scheme runs two alternatives A and B depending on the ratio of η and μ and some small positive input parameters θA and κ.

### Algorithm 1 (Safem(θA,ρB,κ,T0).).

The routine Refine applies the newest vertex bisection [20] and refines the marked triangles to compute the smallest admissible refinement 𝒯+1 of 𝒯 with 𝒯𝒯+1.

The data approximation algorithm approx used in Case B is introduced for separate marking in [8, Section 3.3] or [7, 17, 18] with input tolerance ρBμ2 and values (μ(T):T𝒯).

### 2.4 Axioms and Optimal Convergence

Suppose that there exist universal positive constants Λ1,Λ2,Λ3,Λ4 and 0<ρ2<1 that satisfy (A1)(A4) and (B1)(B2) below. Here and in the following 𝒯^𝕋(𝒯) is an admissible triangulation and refinement of some 𝒯𝕋:=𝕋(𝒯0). Recall the definition of the error estimator in (2.2)–(2.4) and the abbreviation η2(𝒯,):=Mη2(𝒯,M) for all 𝒯.

1. (A1)

Stability: For all 𝒯𝕋 and all 𝒯^𝕋(𝒯),

|η(𝒯^,𝒯𝒯^)-η(𝒯,𝒯𝒯^)|Λ1δ(𝒯,𝒯^).
2. (A2)

Reduction: For all 𝒯𝕋 and all 𝒯^𝕋(𝒯),

η(𝒯^,𝒯^𝒯)ρ2η(𝒯,𝒯𝒯^)+Λ2δ(𝒯,𝒯^).
3. (A3)

Discrete Reliability: For all 𝒯𝕋 and all 𝒯^𝕋(𝒯),

δ2(𝒯,𝒯^)Λ3(η2(𝒯,𝒯𝒯^)+μ2(𝒯)).

Let (𝒯) and (σ) be the output of Safem of Section 2.3 and abbreviate η:=η(𝒯)η(𝒯,𝒯) and η2:=η2(𝒯,𝒯)η2(𝒯):=η(𝒯,𝒯)2, etc.

1. (A4)

Quasi-Orthogonality: For all 0,

k=δ2(𝒯k,𝒯k+1)Λ4σ2.

1. (B1)

Rate s Data Approximation: There exists s>0 such that for Tol>0, 𝒯Tol:=𝚊𝚙𝚙𝚛𝚘𝚡(Tol,μ(T):T𝒯0)𝕋 satisfies

|𝒯Tol|-|𝒯0|Λ5Tol-12sandμ2(𝒯Tol)Tol.
2. (B2)

Quasimonotonicity of μ: For all 𝒯𝕋(𝒯) and all 𝒯^𝕋(𝒯),

μ(𝒯^)μ(𝒯).

### Theorem 1 (Quasi-Optimality [8]).

Suppose (A1)(A4) and (B1)(B2). Then there exists some κ0>0 such that any choice of κ,θA and ρB with

0<κ<κ1:=min{κ0,Λ1-2Λ3-1},0<θA<θ0:=1-κΛ12Λ31+Λ12Λ3,

and 0<ρB<1 implies the following. The output (T)N0 and (σ)N0 of Safem (Algorithm 1) of Section 2.3 satisfies the equivalence

Λ5s+sup0(1+|𝒯|-|𝒯0|)sσΛ5s+supN0(1+N)sminσ(𝕋(N)).

The proof of Theorem 1 is given in [8] and not recalled here. The version of this paper is even slightly simplified in that Λ^3:=0, Λ6:=1, and (𝒯,𝒯^)𝒯𝒯^ are more general in [8] and are not displayed explicitly in this paper.

The data approximation axioms (B1)(B2) are discussed in [8], the results apply verbatim for the setting in this paper. This is exemplified in [8, Section 5] for the mixed FEM and hence not further detailed in this paper.

## 3 Verification of (A1)–(A2)

The analysis of (A1)(A2) follows standard arguments and is outlined here for completeness for the problem at hand with piecewise constant coefficients with little emphasis that the global constants Λ1,Λ2 are bounded by the constant Cjc in the discrete jump control.

## Lemma 1 (Discrete Jump Control [8]).

There exists a universal constant Cjc, which depends on the interior angles in the regular triangulation T and the degree kN0, such that any gPk(T) with jumps

[g]E={(g|T+)|E-(g|T-)|Efor E(Ω) with E=T+T-,g|Efor E(Ω)(K)

satisfies

K𝒯|K|12E(K)[g]EL2(E)2CjcgL2(Ω).

The discrete jump control plus triangle inequalities in Lebesgue spaces and in finite-dimensional Euclidean spaces imply the stability (A1). Throughout this section, let (𝐩RT,uRT) and (𝐩RT^,uRT^) denote the discrete solution with respect to 𝒯𝕋 and its refinement 𝒯^𝕋(𝒯), respectively, and let δ(𝒯,𝒯^) be the distance function (2.4).

## Theorem 2 ((A1) Stability).

Axiom (A1) holds with Λ12:=2(ϱ¯-1+b2)(hmax2+Cjc2) for a global lower bound ϱ¯>0 of the piecewise constant eigenvalues of the coefficient matrix A, the supremum of |b|, the maximal mesh-size hmax, and for the constant Cjc from the discrete jump control of Lemma 1.

## Proof.

The reverse triangle inequality in m for m:=|𝒯𝒯^| over the element contributions implies that

|η(𝒯^,𝒯𝒯^)-η(𝒯,𝒯𝒯^)|2T𝒯𝒯^(η(𝒯^,T)-η(𝒯,T))2.

Each of the terms η(𝒯^,T) and η(𝒯,T) is a norm in 4 of terms, which are Lebesgue norms and so allow for a reverse triangle inequality. This leads to

(η(𝒯^,T)-η(𝒯,T))2|T|gL2(T)2+|T|12E(T)[g]τEL2(E)2

with the abbreviation g:=𝐀-1(𝐩RT^-𝐩RT)+(uRT^-uRT)𝐛P1(𝒯^;2). The sum over all T𝒯𝒯^ involves volume terms and edge jumps. Lemma 1 controls the latter terms and so results in

|η(𝒯^,𝒯𝒯^)-η(𝒯,𝒯𝒯^)|2Cjc2gL2(Ω)2+T𝒯𝒯^|T|gL2(T)2.

The mesh-size is bounded from above and so the right-hand side is bounded by the factor hmax2+Cjc2 times the squared L2 norm of g. A triangle inequality and the bounds on the coefficients show

g22(ϱ¯-1+𝐛2)δ2(𝒯,𝒯^).

## Theorem 3 ((A2) Reduction).

Axiom (A2) holds with ϱ2:=2-14 and Λ2:=Λ1.

## Proof.

For the m refined triangles T𝒯^(K):={T𝒯^:TK} of K𝒯𝒯^, the sum η2(𝒯^,𝒯^𝒯) reads

T𝒯^𝒯(|T|𝐀-1𝐩RT^+uRT^𝐛L2(T)2+|T|12F(T)[𝐀-1𝐩RT^+uRT^𝐛]FτFL2(F)2).

The reverse triangle inequalities in 4m and in Lebesgue spaces over triangles and edges and the abbreviation g from the previous proof plus G:=𝐀-1𝐩RT+uRT𝐛 show

η(𝒯^,𝒯^𝒯)(K𝒯𝒯^T𝒯^(K)(|T|GL2(T)2+|T|12E(T)[G]EτEL2(E)2))12
+(T𝒯^𝒯(|T|gL2(T)2+|T|12E(T)[g]EτEL2(E)2))12.

Since [GτF]F=0 for F^(int(K)) and |T||K|2 for T𝒯^(K), the first term on the right-hand side of the above displayed formula is bounded from above by 2-14η(𝒯,𝒯𝒯^). The remaining term is estimated with Lemma 1 as in the previous proof and so with the same bound on Λ2=Λ1. ∎

## 4 Verification of (A3)

Throughout this section let (𝐩RT,uRT) and (𝐩RT^,uRT^) solve (1.4) and let Π0f and Π^0f denote the L2 orthogonal projection of the right-hand side f onto piecewise constants (P0(𝒯) and P0(𝒯^)) with respect to the triangulation 𝒯 and its refinement 𝒯^, respectively.

The main residual R1 is defined, for any test function 𝐪RT^RT0(𝒯^), by

(4.1)R1(𝐪RT^):=-(𝐀-1𝐩RT+uRT𝐛,𝐪RT^)L2(Ω)+(div𝐪RT^,uRT)L2(Ω).

## Lemma 1.

There exists some qRT^RT0(T^) with norm qRT^H(div,Ω)=1 and

δ(𝒯,𝒯^)R1(𝐪RT^)+Π^0f-Π0f.

## Proof.

The initial mesh-size h0 is sufficiently small throughout this paper to guarantee the existence and stability of the discrete solutions [4, Theorem 4.3]. The stability of the discrete problem (1.4) with respect to the refined triangulation 𝒯^ leads to the existence of (𝐪RT^,vRT^)RT0(𝒯^)×P0(𝒯^) with 𝐪RT^H(div,Ω)+vRT^1 and

δ(𝒯,𝒯^)=(𝐀-1(𝐩RT^-𝐩RT)+(uRT^-uRT)𝐛,𝐪RT^)L2(Ω)-(div𝐪RT^,uRT^-uRT)L2(Ω)
+(div(𝐩RT^-𝐩RT)+γ(uRT^-uRT),vRT^)L2(Ω).

Since (𝐩RT^,uRT^) solves (1.4) with respect to the refined triangulation 𝒯^ and div𝐩RT+γuRT=Π0f (1.4) with respect to 𝒯, this reads

δ(𝒯,𝒯^)=R1(𝐪RT^)+(Π^0f-Π0f,vRT^)L2(Ω).

A Cauchy inequality concludes the proof. ∎

The further analysis of R1 requires a discrete Helmholtz decomposition on a regular triangulation 𝒯 of a (possibly) multi-connected domain Ω. The connectivity components Γ0,Γ1,,ΓJ of Ω are enumerated such that Γ0 denotes the boundary of the unbounded component of 2Ω¯. The modified lowest-order Crouzeix–Raviart space reads

(4.2)CR1(𝒯)={vCR1(𝒯):there exist c1,,cJc0:=0, such that for all j=0,1,,Jand all E(Γj)v(mid(E))=cj}.

(Here and throughout the paper, (Γj) denotes the set of edges on Γj.)

## Lemma 2 (Discrete Helmholtz Decomposition).

For the multi-connected domain Ω the decomposition of piecewise constant vector functions

P0(𝒯;2)=𝐀NCCR1(𝒯)Curl(S1(𝒯))

is orthogonal with respect to the L2 scalar product weighted by A-1 in the sense that

(NCvCR,CurlwC)L2(Ω)=0for all vCRCR1(𝒯) and all wCS1(𝒯).

## Proof.

The discrete Helmholtz decomposition is well known for simply-connected domains and J=0. The general case follows with the same argument by counting triangles and edges |𝒯|+|𝒩|=||+1-J in general; further details are omitted. ∎

The modified Crouzeix–Raviart space is accompanied by a modified Raviart–Thomas space Q(𝒯^) with vanishing integral mean of the normal components over each boundary,

Q(𝒯^):={𝐪RT^RT0(𝒯^):Γj𝐪RT^νds=0 for all j=1,,J}.

The discrete Helmholtz decomposition allows for a characterization of the divergence-free functions.

## Lemma 3 (Discrete Divergence).

The linear operator

div:Q(𝒯^)P0(𝒯^)

is surjective and its kernel is Curl(S1(T^)).

## Proof.

The divergence-free Raviart–Thomas functions in Q(𝒯^) are piecewise constant and allow for a discrete Helmholtz decomposition as in Lemma 2. The decomposition implies that the divergence-free Raviart–Thomas functions in Q(𝒯^) are those in Curl(S1(𝒯^)). Consequently, the kernel of div:Q(𝒯^)P0(𝒯^) has the dimension |𝒩^|-1 for the number |𝒩^| of nodes in 𝒯^. Since the dimension of the vector space Q(𝒯^) is |^|-J for the number |^| of edges in 𝒯^, the range of div:Q(𝒯^)P0(𝒯^) has dimension |^|-|𝒩^|+1-J. Since |𝒯^|+|𝒩^|=|^|+1-J, the range is P0(𝒯^). ∎

One key argument of the reliability analysis is the split of the difference 𝐩RT^-𝐩RT into two parts 𝐩RT^-𝐩RT^ and 𝐩RT^-pRT for some auxiliary solution: Seek (𝐩RT^,uRT^)RT0(𝒯^)×P0(𝒯^) with

(4.3)

(4.3a)Γj𝐩RT^νds=Γj𝐩RTνds,
(4.3b)(𝐀-1𝐩RT^,𝐪RT^)L2(Ω)-(div𝐪RT^,uRT^)L2(Ω)=-(uRT𝐛,𝐪RT^)L2(Ω),
(4.3c)(div𝐩RT^,vRT^)L2(Ω)=(Π0f-γuRT,vRT^)L2(Ω)

hold for all 𝐪RT^Q(𝒯^), for all vRT^P0(𝒯^), and for all j=1,,J.

The solution to (4.3) is recovered from an auxiliary nonconforming problem: Let uCR^CR1(𝒯^) denote the Riesz representation of the functional on the right-hand side of

(𝐀NCuCR^,NCwCR^)L2(Ω)=(Π0f-γuRT,wCR^)L2(Ω)-(uRT𝐛,NCwCR^)L2(Ω)
(4.4)-j=1J(Γj𝐩RTνds)(ΓjwCR^ds)

for all test functions wCR^CR1(𝒯^) in the Hilbert space (CR1(𝒯^),(𝐀NC,NC)L2(Ω)). Here and throughout the paper, denotes the integral mean Γjds:=Γjds/|Γj| for the length |Γj| of the closed polygon Γj.

The equivalence of this with (4.3) implies the unique solvability of (4.3). To verify this, let uCR^ solve (4.4) and, for 𝐱T𝒯^, set

S(T):=T(𝐱-mid(T))𝐀-1(𝐱-mid(T))dx,
(4.5)𝐩RT^(𝐱):=-(𝐀NCuCR^+uRT𝐛)+12(Π0f-γuRT)(𝐱-mid(T)),
(4.6)uRT^(𝐱):=Π0uCR^+S(T)(Π0f-γuRT).

The piecewise constant function S(𝒯) is

S(𝒯)|T=S(T):=14T(𝐱-mid(T))𝐀-1((𝐱-mid(T))dx

in T𝒯.

## Lemma 4.

The pair (pRT^,uRT^) from (4.5)–(4.6) defines the unique solution to (4.3).

## Proof.

The proof imitates that of [4, Theorem 4.2] and generalizes it to multiply connected domains. The arguments therein confirm the continuity of normal components along the interior edges and prove 𝐩RT^RT0(𝒯^)H(div,Ω). The present situation involves the connectivity components Γ1,,ΓJ of the boundary Ω and requires a little modification in the proof that (𝐩RT^,uRT^) indeed satisfies (4.3). Set

cj:=ΓjuCR^ds

and substitute 𝐩RT^ from (4.5) before a piecewise integration by parts shows, for any test function 𝐪RT^Q(𝒯^), that

(𝐩RT^+uRT𝐛,𝐪RT^)L2(Ω)-(div𝐪RT^,Π0uCR^)L2(Ω)
=(12𝐀-1(Π0f-γuRT)(-mid(𝒯)),𝐪RT^)L2(Ω)-j=1JcjΓj(𝐪RT^ν)ds.

The first term on the right-hand side already appears in [4, p. 567, lines 1–2] and is rewritten as

2((Π0f-γuRT)𝐀-1(-mid(𝒯)),𝐪RT^)L2(Ω)=(Π0f-γuRT,S(𝒯)div𝐪RT^)L2(Ω).

This term combines with -(div𝐪RT^,Π0uCR^)L2(Ω) in the aforementioned equality and leads with (4.6) to -(div𝐪RT^,uRT^)L2(Ω). Recall that c0=0 and that Γj𝐪RT^νds=0 for j=1,,J because of 𝐪RT^Q(𝒯^). This calculation leads to equation (4.3b). Since 𝐩RT^RT0(𝒯^)H(div,Ω), the definition of (4.5) immediately proves (4.3c). Rewrite (4.5) to obtain an identity for 𝐀NCuCR^ and utilize this in (4.4). This leads to an identity, which allows a piecewise integration by parts and then results in

j=1JΓjwCR^(𝐩RT^-𝐩RT)νds=0for all wCR^CR1(𝒯^).

The design of wCR^CR1(𝒯^) with piecewise integral means along the boundary edges of Γj which are constant for each j=1,,J proves (4.3a).

This concludes the proof of the existence of a discrete solution to (4.3a)–(4.3c) and it remains to show the uniqueness of a discrete solution. This follows from the trivial solution to the homogeneous system

(𝐀-1𝐩RT^,𝐪RT^)L2(Ω)-(div𝐪RT^,uRT^)L2(Ω)=0=(div𝐩RT^,vRT^)L2(Ω)

for all 𝐪RT^Q(𝒯^), vRT^P0(𝒯^) plus Γj𝐩RT^νds=0 for j=0,,J. Given an arbitrary solution (𝐩RT^,uRT^) to this discrete homogeneous problem, let 𝐪RT^=𝐩RT^ to find 𝐩RT^=0. It follows (div𝐪RT^,uRT^)L2(Ω)=0 for all test functions. Lemma 3 allows for a test function 𝐪RT^ with div𝐪RT^=uRT^ and so uRT^=0. This concludes the proof of the uniqueness of the solution of the homogeneous system. ∎

## Lemma 5.

The test function qRT^RT0(T^) in Lemma 1 can be selected additionally to satisfy Π0divqRT^=0 and ΓjqRT^νds=0 for all j=0,1,,J.

## Proof.

The first equation in (1.4) shows that RT0(𝒯) belongs to the kernel of R1 from (4.1). Hence we may and will replace test function 𝐪RT^RT0(𝒯^) in Lemma 1 by 𝐪RT^-𝐪RT for some appropriate 𝐪RTRT0(𝒯). The naive choice of the Fortin interpolation 𝐪RT of 𝐪RT^ leads to the additional properties but leaves open the subtle question of the uniform bound 𝐪RTH(div,Ω) in terms of 𝐪RT^H(div,Ω)1. This proof utilizes the MFEM solution (𝐪RT,vRT)RT0(𝒯)×P0(𝒯) to the modified Poisson model problem with right-hand side Π0div𝐪RT^,

(4.7)

(4.7a)ΓJ𝐪RTνds=Γj𝐪RT^νdsfor all j=1,,J,
(4.7b)(𝐪RT,𝐳RT)L2(Ω)=(div𝐳RT,vRT)L2(Ω)for all 𝐳RTRT0(𝒯),
(4.7c)(div𝐪RT,wRT)L2(Ω)=(Π0div𝐪RT^,wRT)L2(Ω)for all wRTP0(𝒯).

Equations (4.7) are a particular version of (4.3) and the equivalence of Lemma 4 applies here as well. This implies the unique existence of (𝐪RT,vRT)RT(𝒯)×P0(𝒯) and leads the desired bound

𝐪RTH(div,Ω)𝐪RT^H(div,Ω)1.

Recall that (𝐩RT,uRT) and (𝐩RT^,uRT^) solve (1.4) with respect to the triangulation 𝒯 and its refinement 𝒯^, respectively, and that (𝐩RT^,uRT^) solves (4.3). Define the modified residual R2 by

(4.8)R2(𝐪RT^):=-(𝐀-1𝐩RT^+uRT𝐛,𝐪RT^)L2(Ω)for any test function 𝐪RT^RT0(𝒯^).

Lemma 5 and the modified test-function 𝐪RT^ show

δ(𝒯,𝒯^)-Π^0f-Π0fR1(𝐪RT^)=-(𝐀-1𝐩RT+uRT𝐛,𝐪RT^)L2(Ω)
(4.9)R2(𝐪RT^)+𝐩RT^-𝐩RTH(div,Ω).

The divergence-free term 𝐩RT^-𝐩RT is controlled in the subsequent lemma.

## Lemma 6.

It holds pRT^-pRTH(div,Ω)η(T,TT^).

## Proof.

Since div𝐩RT^=Π0f-γuRT=div𝐩RT implies div(𝐩RT^-𝐩RT)=0, it follows that 𝐩RT^-𝐩RT is piecewise constant and its discrete Helmholtz decomposition leads to αCR^CR1(𝒯^) and βC^S1(𝒯^) with

𝐩RT^-𝐩RT=𝐀NCαCR^+CurlβC^.

This and the L2 orthogonality (𝐩RT^-𝐩RT)NCαCR^ show that αCR^0. Consequently,

(4.10)𝐩RT^-𝐩RTH(div,Ω)2=(𝐩RT^-𝐩RT,𝐀-1CurlβC^)L2(Ω).

Given any node z𝒩 in the coarse triangulation 𝒯, the Scott–Zhang quasi-interpolation [19] defines βC(z) by a selection of an edge E(z) with vertex z and evaluates some weighted integral of βC^ along E(z). Select the edge E(z)^ if possible to obtain a Scott–Zhang quasi-interpolation βCS1(𝒯) of βC^ with βC^-βC=0 a.e. in (𝒯𝒯^) plus the local approximation and stability properties. For any edge E of length hE and its neighborhood Ω(E):=z𝒩(E)ωz for the nodal patches ωz, the latter properties and discrete trace inequalities result in

(4.11)βC^-βCL2(E)hE12βC^H1(Ω(E)).

The weak formulation (1.4) with 𝐪RT=CurlβCRT0(𝒯) and equation (4.3b) with 𝐪RT^:=CurlβCRT0(𝒯^) show (𝐩RT^-𝐩RT,𝐀-1CurlβC)L2(Ω)=0. Hence (4.10) is (𝐩RT^-𝐩RT,𝐀-1Curl(βC^-βC))L2(Ω). The test function 𝐪RT^:=Curl(βC^-βC) in equation (4.3b) shows that (4.10) is (𝐀-1𝐩RT+uRT𝐛,Curl(βC-βC^))L2(Ω). This and a piecewise integration by parts leads to

𝐩RT^-𝐩RTH(div,Ω)2=E^E(βC-βC^)[𝐀-1𝐩RT+uRT𝐛]EτEds-T𝒯^T(βC-βC^)curl(𝐀-1𝐩RT+uRT𝐛)ds.

The piecewise curl of the low-order Raviart–Thomas finite element functions vanishes and so do all summands in the last term. Since βC^-βC=0 along any edge E^, this proves

𝐩RT^-𝐩RTH(div,Ω)2E^[𝐀-1𝐩RT+uRT𝐛]EτEL2(E)βC^-βCL2(E).

The combination with estimate (4.11) for βC^-βCL2(E) and the bound βC^=CurlβC^𝐩RT^-𝐩RTH(div,Ω) imply

(4.12)𝐩RT^-𝐩RTH(div,Ω)2E^hE[𝐀-1𝐩RT+uRT𝐛]EτEL2(E)2.

A rearrangement with the triangle-oriented error estimator concludes the proof. ∎

## Lemma 7.

The test function qRT^ from Lemma 5 and the residual R2 from (4.8) satisfy

R2(𝐪RT^)η(𝒯,𝒯𝒯^)+𝐩RT^-𝐩RTH(div,Ω).

## Proof.

Let -mid(𝒯^)P1(𝒯^;2) abbreviate the function x-mid(T) for xT𝒯^ and consider the test function 𝐪RT^=Π^0𝐪RT^+12div𝐪RT^(-mid(𝒯^))P1(𝒯^;2) from Lemma 5. The piecewise constant part Π^0𝐪RT^ allows a discrete Helmholtz decomposition

Π^0𝐪RT^=𝐀NCvCR^+CurlβC^

from Lemma 2 for some vCR^CR1(𝒯^) and some βC^S1(𝒯^). Altogether, there are three contributions of

R2(𝐪RT^)=R2(𝐀NCvCR^)+R2(CurlβC^)+R2(12div𝐪RT^(-mid(𝒯^))).

(i) The representation of 𝐩RT^ in (4.5) is utilized in the first contribution and leads to

R2(𝐀NCvCR^)=-(𝐩RT^+uRT𝐛,NCvCR^)L2(Ω)=(NCuCR^,𝐀NCvCR^)L2(Ω).

On the other hand, 𝐀NCvCR^=Π^0𝐪RT^-CurlβC^ and CurlβC^ is L2 orthogonal to NCuCR^. This and an integration by parts prove

R2(𝐀NCvCR^)=(NCuCR^,𝐪RT^)L2(Ω)=-(uCR^,div𝐪RT^)L2(Ω).

Recall that div𝐪RT^ is L2 orthogonal onto P0(𝒯) and so

R2(𝐀NCvCR^)=((Π0-Π^0)uCR^,div𝐪RT^)L2(Ω)(Π^0-Π0)uCR^div𝐪RT^.

Note that (Π^0-Π0)uCR^ vanishes a.e. in T𝒯𝒯^ and a piecewise discrete Poincaré inequality shows (Π^0-Π0)uCR^h𝒯NCuCR^L2(Ω) with the interior Ω of the domain (𝒯𝒯^). Consequently,

R2(𝐀NCvCR^)h𝒯NCuCR^L2(Ω).

(ii) Observe from (4.3b) and CurlβC^Q(𝒯^) that R2(CurlβC^)=0.

(iii) Since div𝐪RT^ vanishes outside the set Ω (for div𝐪RT^=Π0div𝐪RT^ on T𝒯𝒯^) and the weight satisfies |-mid(𝒯^)|h𝒯, the term R2(12div𝐪RT^(-mid(𝒯^))) is equal to

-12(𝐀-1𝐩RT^+uRT𝐛,(-mid(𝒯^))div𝐪RT^)L2(Ω)h𝒯(𝐀-1𝐩RT^+uRT𝐛)L2(Ω).

In conclusion of (i)–(iii), it follows that

(4.13)R2(𝐪RT^)h𝒯NCuCR^L2(Ω)+h𝒯(𝐀-1𝐩RT^+uRT𝐛)L2(Ω).

On the other hand, the representation formula (4.5) shows that

h𝒯(𝐀-1𝐩RT^+uRT𝐛-NCuCR^)L2(Ω)=h𝒯(Π0f-γuRT)𝐀-1(-mid(𝒯))L2(Ω).

With a lower bound ϱ¯ of the smallest eigenvalue of A and with div𝐩RT=Π0f-γuRT from (1.4),

h𝒯(𝐀-1𝐩RT^+uRT𝐛-NCuCR^)L2(Ω)ϱ¯-1h𝒯2div𝐩RTL2(Ω).

For each T𝒯, hTdiv𝐩RTL2(T)(1-Π0)𝐩RTL2(T) and an inverse estimate plus a uniform upper bound ϱ¯ of the eigenvalues of 𝐀 lead to

ϱ¯h𝒯2div𝐩RTL2(T)ϱ¯h𝒯(1-Π0)(𝐀-1𝐩RT)L2(T)
h𝒯(1-Π0)(𝐀-1𝐩RT)L2(T)h𝒯(𝐀-1𝐩RT+uRT𝐛)L2(T).

The combination of the previously displayed estimates results in

h𝒯(𝐀-1𝐩RT^+uRT𝐛-NCuCR^)L2(Ω)h𝒯(𝐀-1𝐩RT+uRT𝐛)L2(Ω).

This and (4.13) plus some triangle inequalities imply

R2(𝐪RT^)h𝒯(𝐀-1𝐩RT+uRT𝐛)L2(Ω)+h𝒯(𝐀-1𝐩RT^+uRT𝐛)L2(Ω)
h𝒯(𝐩RT^-𝐩RT)L2(Ω)+h𝒯(𝐀-1𝐩RT+uRT𝐛)L2(Ω).

Since h𝒯h01, this concludes the proof. ∎

## Theorem 8 ((A3) Discrete Reliability).

Under the overall assumption that h0 is sufficiently small, there exists some universal constant Λ3, which depends on the global lower and upper bounds of the eigenvalues of A and on the universal stability constant of the discrete problems and on the shape-regularity in T such that the following holds. The discrete solutions (pRT,uRT) and (pRT^,uRT^) of (1.4) with respect to the triangulation T and its refinement T^ satisfy

Λ3-1δ2(𝒯,𝒯^)η2(𝒯,𝒯𝒯^)+μ2(𝒯)-μ2(𝒯^).

## Proof.

This follows from the combination of Lemma 5, 6, and 7 with (4.9). ∎

Axiom (A3) implies the reliability result of [4, Theorem 5.5].

## Corollary 9 (Reliability).

The solution (p,u) to (1.3) and the solution (pRT,uRT) to (1.4) satisfy

(𝐩-𝐩RT,u-uRT)H(div,Ω)×L2(Ω)Λ312σ(𝒯).

## Proof.

This follows from (A3) for a fixed triangulation 𝒯 and a sequence of its successive uniform refinements 𝒯^ as then the maximal mesh-size in 𝒯^ tends to zero and standard estimates show convergence of (𝐩RT^,uRT^) to (𝐩,u) in the norm of H(div,Ω)×L2(Ω). ∎

## 5 Verification of (A4)

The following lemma proves the supercloseness property (1.5) of Π0u to the mixed solution uRT with a duality argument. For given right-hand side gL2(Ω), the dual problem seeks ϕH01(Ω) with

(5.1)ϕ:=-div(𝐀ϕ)+𝐛ϕ+γϕ=g.

Under the overall assumption that is injective, it follows that and its dual are isomorphisms between H01(Ω) and H-1(Ω).

The reduced elliptic regularity of the leading elliptic part -div(𝐀) leads to higher regularity, that is, there exist α with 0<α1 and Creg< with

(5.2)|ϕ|H1+α(Ω)Cregg

for any right-hand side gL2(Ω) with the solution ϕ to (5.1) (see [12, Sections 5 and 14]).

The supercloseness (1.5) is discussed in the introduction and (unlike the remaining results of this paper) holds without any assumption on the initial mesh-size as long as (1.1) is injective and (1.4) has a solution.

## Lemma 1 (Supercloseness).

The solution (p,u) to (1.3) and the solution (pRT,uRT) to (1.4) satisfy (1.5).

## Proof.

The dual problem (5.1) and its solution ϕH01(Ω)H1+α(Ω) for the right-hand side g=Π0u-uRTP0(𝒯) lead to 𝐪:=𝐀ϕH(div,Ω)Lt(Ω;2) for some t>2. This allows the application of the Fortin interpolation operator IF [3, pp. 107–109] with the commutative property Π0div𝐪=divIF𝐪. This and ϕh:=Π0ϕ result in

g2=-(g,div𝐪)L2(Ω)+(g,𝐛ϕ+γϕ)L2(Ω)
=-(u-uRT,divIF𝐪)L2(Ω)+(g,Π0(𝐛ϕ)+γϕh)L2(Ω).

Equations (1.3)–(1.4) show Π0div(𝐩-𝐩RT)=-γg and

(u-uRT,divIF𝐪)L2(Ω)=(𝐀-1(𝐩-𝐩RT)+(u-uRT)𝐛,IF𝐪)L2(Ω).

The combination of the aforementioned identities in the first step and the identity ϕ=𝐀-1𝐪 in the second step plus an integration by parts prove that

(5.3)g2=(𝐀-1(𝐩-𝐩RT)+(u-uRT)𝐛,𝐪-IF𝐪)L2(Ω)+(div(𝐩-𝐩RT),ϕ-ϕh)L2(Ω)-(u-uRT,(1-Π0)(𝐛ϕ))L2(Ω).

The error estimates of the Fortin interpolation plus piecewise Poincaré inequalities and the reduced elliptic regularity (5.2) of the dual problem (5.1) imply

𝐪-IF𝐪hmaxα|𝐪|Hα(Ω)hmaxαg,
ϕ-Π0(ϕ)hmaxαϕH1+α(Ω)hmaxαg,
ϕ-ϕhhmaxϕhmaxg.

The application of these approximation properties to (5.3) concludes the proof. ∎

Recall that α>0 is the positive extra regularity parameter, which exclusively depends on the domain and on the coefficients, and the maximal initial mesh-size h0 is the maximal mesh-size in 𝒯0 (whence in all 𝒯𝕋) and Λ3 is from (A3).

## Theorem 2 ((A4ϵ) Quasi-Orthogonality).

There exists a constant Λ4< such that for sufficiently small h0, any ,mN0 satisfy

Λ3-1k=+mδ2(𝒯k,𝒯k+1)2σ2+Λ4h02αk=+1+mσ2.

## Proof.

Recall 𝐩-𝐩H(div,Ω)2:=𝐀-12(𝐩-𝐩)2+div(𝐩-𝐩)2 and abbreviate

e2:=𝐩-𝐩H(div,Ω)2+u-u2andδ,+12:=𝐩+1-𝐩H(div,Ω)2+u+1-u2.

Elementary algebra plus (1.3)–(1.4) with respect to the triangulation 𝒯+1, each with the test function (𝐩+1-𝐩,u+1-u), eventually shows that δ,+12+e+12-e2 is equal to

(5.4)2(u-u+1,𝐛(𝐩+1-𝐩)+(γ-1)div(𝐩+1-𝐩)-(u+1-u))L2(Ω).

The factor u-u+1 in this L2 scalar product is split into Π+1u-u+1 and u-Π+1u with the L2 projection Π+1 onto P0(𝒯+1). Lemma 1 applies on the level of 𝒯+1 and controls the L2 norm

Π+1u-u+1Ch0αe+1.

This and the Cauchy inequality control the first contribution in (5.4) by h0αe+1δ,+1 times the constant 2C(𝐀-12𝐛2+γ-12+1)12. The remaining second contribution in term (5.4) is the L2 scalar product of 2(u-Π+1u)𝐛 with (1-Π+1)(𝐩+1-𝐩) (recall that 𝐛=𝐀-1𝐛 is a constant vector with length |𝐛|=𝐛). The Raviart–Thomas function 𝐩+1-𝐩RT0(𝒯+1) allows in 2D for

(1-Π+1)(𝐩+1-𝐩)(x)=12(x-mid(T))div(𝐩+1-𝐩)

at xT𝒯+1 and so

2((u-Π+1u)𝐛,𝐩+1-𝐩)L2(Ω)23|𝐛|h0u-u+1div(𝐩+1-𝐩)23|𝐛|h0e+1δ,+1.

The combination of the two estimates for the two contributions in (5.4) leads (with h0h0α for α1 and h01) to

Λ4:=(23𝐛+2C(𝐀-12𝐛2+γ-12+1)12)2,
δ,+12+e+12-e2h0α(Λ4)12e+1δ,+112δ,+12+12Λ4h02αe+12.

Provided that Λ4h02α2, the sum of the above inequalities over different levels shows

k=+mδk,k+122e2+Λ4h02αk=+1+mek2.

This and Corollary 9 with ek2Λ3σk2 conclude the proof. ∎

## Corollary 3 ((A4) Quasi-Orthogonality).

Axiom (A4) holds for sufficiently small initial mesh-sizes h0.

## Proof.

This follows from [8, Theorem 3.1] and (A4ϵ) for a sufficiently small ϵ:=Λ3Λ4h02α<1-ρ12Λ12 with parameters 0<ρ12<1 and Λ12>0 of the reduction result (A12) in [8, Theorem 4.1]. ∎

## Remark 4 (Generalizations).

Several arguments in this section apply to other mixed finite element methods as well but the second contribution in (5.4) solely applies to the Raviart–Thomas mixed finite element family. The restriction on the smallness of h0 in Theorem 2 can be circumvented by the quasi-monotonicity (cf. [8] for the concept) but is required in Corollary 3.

## 6 Numerical Experiments

This section is devoted to numerical experiments to investigate the influence of the critical parameters h0, θA, and κ and the practical performance of the adaptive algorithm Safem. After a few remarks on the implementation, three examples on the L-shaped domain are displayed with smooth or discontinuous right-hand sides, before some overall observations conclude the paper.

### 6.1 Numerical Realization

The data approximation is realized by the Thresholding Second Algorithm (TSA) of [1] followed by the closure algorithm to output a shape-regular triangulation.

### Algorithm 2 (TSA.).

The realization from [17, 18] is slightly modified in the Approx algorithm of [1] through a parameter ϑ=1-10-6<1 in the computation of . The functional μ~(T) in TSA is a weighted error functional, which depend on the values of μ(T) and μ~(T) on the parent triangle TP as well as on the siblings of TP, cf. [1, 17, 18] for more details and the explicit formulas.

The non-homogeneous boundary data in Section 6.2 are not met in the theoretical part of this paper, which is simplified to homogeneous boundary conditions. The first example with known solution requires inhomogeneous boundary data on Ω with the modified jump-term (𝐀-1𝐩RT+uRT𝐛)|EτE+us|E along the boundary edge EΩ with the prescribed boundary values u and its tangential derivative us on E.

#### Short Notation.

The abbreviation error ε (resp. estimator σ) refers to the left-(resp. right-)hand side of (2.1).

### 6.2 Continuous Right-Hand Side with Known Corner Singularity

The coefficients 𝐀=I, 𝐛=(1,1) and γ=-2 on the L-shaped domain Ω=(-1,1)2[0,1)×(-1,0] allow in (1.1) for a right-hand side f (computed by (1.1) for given u) and inhomogeneous Dirichlet boundary data (taken from u) for

u(r,θ)=r23sin(2θ3)

in polar coordinates (r,θ) centered at the origin. The initial mesh 𝒯0 consists of 24 congruent right-isosceles triangles from a criss refinement of the three sub-squares and 28 degrees of freedom.

### Figure 1

Convergence history plot for different values of θA in Section 6.2.

Figure 1 displays the outcome of Safem (Algorithm 1) with κ=1, ρB=12, and various values of θA=0.1, 0.3, and 0.5. Despite the condition of Theorem 1 on the smallness of the bulk parameter θA, all displayed values result in an improved optimal empirical convergence rate 12 compared to uniform mesh-refinement with a known suboptimal convergence rate 13. The initial mesh is relatively coarse but all convergence rates are visible right from the beginning, the pre-asymptotic range is not visible. This examples allows the computation of the errors and the estimators and the equivalence is visible throughout with a the expected behavior. In the displayed experiment with κ=1 only Case A of Safem applies as the right-hand side f is continuous. In case of κ0.1 for instance, only Case B runs in Safem for a very long computational range.

### 6.3 Constant Right-Hand Side

The coefficients 𝐀=0.1I, 𝐛=(1,2) and γ=-4 on the L-shaped domain Ω with constant right-hand side f1 lead in (1.1) to an unknown weak solution uH01(Ω). Figure 2 (left) displays the output 𝒯13 of through Safem with κ=1 and ρB=12=θA. Besides the local mesh-refining at the re-entering corner, some layer of refinement are visible along some part of the boundary, that mimics a singular perturbed situation with 𝐀=ϵI for very small ϵ. Undisplayed numerical experiments for ϵ=0.01 have confirmed this observation even stronger.

### Figure 2

Triangulation 𝒯13 on level 13 in Section 6.3 with 𝚗𝚍𝚘𝚏=2254 (left) and Section 6.4 with 𝚗𝚍𝚘𝚏=4485 and ϵ=1 (right).

### Figure 3

Convergence history plot for different initial input triangulations in Safem in Section 6.3ϵ=0.1.

The convergence history plot of Figure 3 displays the estimators σ as functions of the number of degrees of freedom for various initial meshes, namely for 𝒯0 as described in the previous subsection and also for an initial mesh red(𝒯0) (of 𝒯0 from the previous subsection) with 𝚗𝚍𝚘𝚏=140; red-refinement means the division of each triangle into four congruent sub-triangles by connecting its edges’ midpoints with straight lines. This is plotted under the label σ(uniform) and shows the expected suboptimal empirical convergence rate. Those red-refined triangulations, e.g., red2(𝒯0) (with 𝚗𝚍𝚘𝚏=544) with two red-refinements and red3(𝒯0) (with 𝚗𝚍𝚘𝚏=2240) for three, serve as initial triangulations in the input of Safem and Figure 3 displays the respective convergence history plots.

The numerical experiment with the coarsest initial mesh 𝒯0 displays a pre-asymptotic range up to 1000 degrees of freedom. The finer initial triangulations lead to a much smaller pre-asymptotic range with a rapid decrease through a strong local mesh-refinement until the convergence rate of the other adaptive mesh-refinements is met. For the displayed parameter κ=1, solely the Case A runs in Safem.

The undisplayed numerical experiment for a smaller parameter ϵ=0.01 leads to a much larger pre-asymptotic domain with a systematic error reduction only for a fine initial mesh red3(𝒯0).

### 6.4 Piecewise Constant Right-Hand Side

Given the constant coefficients with 𝐀=ϵI and the domain as in the previous subsection, the right-hand side f for this example is piecewise constant with the values ±1 and the value -1 exactly on the two squares ω:=(25,35)2(35,45)2 not aligned to the triangulations; f|Ωω1 and f|ω-1. Figure 2 (right) displays the output 𝒯13 of Safem with κ=1=ϵ and ρB=12=θA with two squares ω visible by local mesh-refinements along ω to resolve the discontinuity of the right-hand side f (recall that f is discontinuous at triangles that intersect ω). Cases A and B alternate in Safem for this example with ϵ=1 for the resolution of the discontinuities of the right-hand side at ω.

The convergence history plot for this example is not displayed as it looks very similar to that of the previous subsection (namely Figure 3) although the reasons for a larger pre-asymptotic range might be different.

For smaller values of ϵ, the triangulations look more like the picture in Figure 2 (left) from the previous subsection and solely Case A runs in Safem. Even for ϵ=0.1 (and more so for ϵ=0.01), the mesh-refinement displays boundary layers and no longer the discontinuities along ω.

### 6.5 Conclusions

The overall impression from the displayed and undisplayed numerical experiments is that the algorithm Safem is very robust such that the choice of θA, ρB, κ in the asymptotic convergence regime with an observed optimal convergence rate: The values θA=12=ρB and κ=1 can be recommended throughout the examples of this paper. The condition on a sufficiently fine initial mesh 0<h01 dramatically influences the pre-asymptotic behavior. Although the examples in Subsection 6.3 and 6.4 are very different in the right-hand side, the stability of the discrete system is identical. The finer the initial mesh, the smaller is the pre-asymptotic range in particular for 𝐀=ϵI with very small ϵ (e.g., ϵ=0.01). This paper exploits the situation when solely is injective and then 0<h01 appears necessary for the well-posedness of the discrete systems and has to be monitored in practise. It is conjectured that this dominates the difficulty of choosing an appropriate initial triangulation in Safem.

Funding source:

Funding statement: The research has been supported by the Deutsche Forschungsgemeinschaft in the Priority Program 1748 under the project Foundation and application of generalized mixed FEM towards nonlinear problems in solid mechanics (CA 151/22-1). The research was carried out while the second author (Asha K. Dond) enjoyed the support of the Priority Program 1748 Reliable simulation techniques in solid mechanics. Development of non-standard discretization methods, mechanical and mathematical analysis for a research visit at the Humboldt-Universität zu Berlin, Germany.

# Acknowledgements

The work has been written, while the first author enjoyed the hospitality of the Hausdorff Research Institute of Mathematics in Bonn, Germany, during the Hausdorff Trimester Program Multiscale Problems: Algorithms, Numerical Analysis and Computation. The authors thank Rui Ma for a careful reading of the manuscript and her valuable remarks.

### References

[1] P. Binev, W. Dahmen and R. DeVore, Adaptive finite element methods with convergence rates, Numer. Math. 97 (2004), no. 2, 219–268. Search in Google Scholar

[2] P. Binev and R. DeVore, Fast computation in adaptive tree approximation, Numer. Math. 97 (2004), no. 2, 193–217. Search in Google Scholar

[3] D. Boffi, F. Brezzi and M. Fortin, Mixed Finite Element Methods and Applications, Springer Ser. Comput. Math. 44, Springer, Heidelberg, 2013. Search in Google Scholar

[4] C. Carstensen, A. K. Dond, N. Nataraj and A. K. Pani, Error analysis of nonconforming and mixed FEMs for second-order linear non-selfadjoint and indefinite elliptic problems, Numer. Math. 133 (2016), no. 3, 557–597. Search in Google Scholar

[5] C. Carstensen, M. Feischl, M. Page and D. Praetorius, Axioms of adaptivity, Comput. Math. Appl. 67 (2014), no. 6, 1195–1253. Search in Google Scholar

[6] C. Carstensen and R. H. W. Hoppe, Error reduction and convergence for an adaptive mixed finite element method, Math. Comp. 75 (2006), no. 255, 1033–1042. Search in Google Scholar

[7] C. Carstensen and H. Rabus, An optimal adaptive mixed finite element method, Math. Comp. 80 (2011), no. 274, 649–667. Search in Google Scholar

[8] C. Carstensen and H. Rabus, Axioms of adaptivity with separate marking for data resolution, SIAM J. Numer. Anal. 55 (2017), no. 6, 2644–2665. Search in Google Scholar

[9] J. M. Cascon, C. Kreuzer, R. H. Nochetto and K. G. Siebert, Quasi-optimal convergence rate for an adaptive finite element method, SIAM J. Numer. Anal. 46 (2008), no. 5, 2524–2550. Search in Google Scholar

[10] H. Chen, X. Xu and R. H. W. Hoppe, Convergence and quasi-optimality of adaptive nonconforming finite element methods for some nonsymmetric and indefinite problems, Numer. Math. 116 (2010), no. 3, 383–419. Search in Google Scholar

[11] L. Chen, M. Holst and J. Xu, Convergence and optimality of adaptive mixed finite element methods, Math. Comp. 78 (2009), no. 265, 35–53. Search in Google Scholar

[12] M. Dauge, Elliptic Boundary Value Problems on Corner Domains. Smoothness and Asymptotics of Solutions, Lecture Notes in Math. 1341, Springer, Berlin, 1988. Search in Google Scholar

[13] A. K. Dond, N. Nataraj and A. K. Pani, Convergence of an adaptive lowest-order Raviart–Thomas element method for general second-order linear elliptic problems, IMA J. Numer. Anal. 37 (2017), no. 2, 832–860. Search in Google Scholar

[14] W. Dörfler, A convergent adaptive algorithm for Poisson’s equation, SIAM J. Numer. Anal. 33 (1996), no. 3, 1106–1124. Search in Google Scholar

[15] M. Feischl, T. Führer and D. Praetorius, Adaptive FEM with optimal convergence rates for a certain class of nonsymmetric and possibly nonlinear problems, SIAM J. Numer. Anal. 52 (2014), no. 2, 601–625. Search in Google Scholar

[16] K. Mekchay and R. H. Nochetto, Convergence of adaptive finite element methods for general second order linear elliptic PDEs, SIAM J. Numer. Anal. 43 (2005), no. 5, 1803–1827. Search in Google Scholar

[17] H. Rabus, Quasi-optimal convergence of AFEM based on separate marking, Part I, J. Numer. Math. 23 (2015), no. 2, 137–156. Search in Google Scholar

[18] H. Rabus, Quasi-optimal convergence of AFEM based on separate marking, Part II, J. Numer. Math. 23 (2015), no. 2, 157–174. Search in Google Scholar

[19] L. R. Scott and S. Zhang, Finite element interpolation of nonsmooth functions satisfying boundary conditions, Math. Comp. 54 (1990), no. 190, 483–493. Search in Google Scholar

[20] R. Stevenson, Optimality of a standard adaptive finite element method, Found. Comput. Math. 7 (2007), no. 2, 245–269. Search in Google Scholar

[21] R. Stevenson, The completion of locally refined simplicial partitions created by bisection, Math. Comp. 77 (2008), no. 261, 227–241. Search in Google Scholar